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SUMMARY 


LEWICE is an ice accretion prediction code that applies a time-stepping 
procedure to calculate the shape of an ice accretion. The potential flow field 
is calculated in LEWICE using the Douglas Hess-Smith 2-D panel code 
(S24Y). This potential flow field is then used to calculate the trajectories 
of particles and the impingment points on the body. These calculations 
are performed to determine the distribution of liquid water impinging on 
the body, which then serves as input to the icing thermodynamic code. 
The icing thermodynamic model is based on the work of Messinger, but 
contains several major modifications and improvements. This model is 
used to calculate the ice growth rate at each point on the surface of the 
geometry. By specifying an icing time increment, the ice growth rate can 
be interpreted as an ice thickness which is added to the body, resulting in 
the generation of new coordinates. This procedure is repeated, beginning 
with the potential flow calculations, until the desired icing time is reached. 

The operation of LEWICE is illustrated through the use of five ex- 
amples. These examples are representative of the types of applications 
expected for LEWICE. All input and output is discussed, along with many 
of the diagnostic messages contained in the code. Several error conditions 
that may occur in the code for certain icing conditions are identified, and 
a course of action is recommended. 

LEWICE has been used to calculate a variety of ice shapes, but should 
still be considered a research code. The code should be exercised further to 
identify any shortcomings and inadequacies. Any modifications identified 
as a result of these cases, or of additional experimental results, should be 
incorporated into the model. Using it as a test bed for improvements to 
the ice accretion model is one important application of LEWICE. 


Chapter 1 

BACKGROUND 

The evaluation of both commercial and military flight systems in icing 
conditions has become important in the design and certification phases 
of system development. These systems have been evaluated in flight in 
natural icing, in a simulated cloud produced by a leading aircraft, and in 
ground test facilities. All icing testing is relatively expensive, and each test 
technique, i.e., flight or ground testing, has operational limitations which 
limit the range of icing conditions that can be evaluated. It would benefit 
the aircraft or flight system manufacturer to be able to analytically predict 
the performance of the system for a range of icing conditions. 

This first step in the prediction of the performance characteristics is 
the determination of the location, size, and shape of the ice that will form. 
An analytical ice accretion model would allow the evaluation of a wide 
range of proposed test conditions to identify those that will be most critical 
to the flight system. This could substantially reduce the amount of test 
time required to adequately evaluate a system and increase the quality 
and confidence level of the final evaluation. The analytically predicted 
ice accretion could also serve as the input to an advanced aerodynamic or 
system performance code to allow more complete evaluation in the design 
phases of the system. For these reasons, several analytical ice accretion 
prediction methods have been developed by various investigators. 

The most well-known are those of Ackley and Templeton 1 , Lowzowski 2 , 
Hankey and Kirchner 3 , and Cansdale and Gent 4 . All apply essentially the 
same physical model of the ice accretion process but differ in the manner 
that the surface properties, for example, the local collection efficiency and 
convective heat transfer coefficient, are calculated and allowed to vary over 
the surface. Also, severed of these models were restricted to the simulation 
of the icing of a cylinder. One of the major inadequacies of these models 
is that they do not account for any of the time-dependent aspects of the 
ice accretion process. In general, the ice accretion rate is calculated as 
a function of position on the airfoil and projected at a constant rate to 
approximate a finite growth over a prescribed period of time. Therefore, 
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any sensitivity of the ice accretion process to the changing iced airfoil shape 
is not included. 

The purpose of the current study was to develop a time-dependent, an- 
alytical model of the ice accretion process that could be used to predict the 
shape of the ice accretion that would form on an arbitrary two-dimensional 
geometry when exposed to icing conditions. The development of the com- 
puter code (LEWICE) was begun by the University of Dayton Research 
Institute (UDRI) under contract to NASA Lewis Research Center. The 
results of this study are described in Reference 5. Development of the code 
was then continued at the NASA Lewis Research Center under NASA fund- 
ing until October 1984 when funding by the Federal Aviation Administra- 
tion (FAA) was begun. This document describes the results of the current 
development effort and contains the information necessary to apply the ice 
accretion prediction method to practical icing problems. 
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Chapter 2 
INTRODUCTION 

The computer code, LEWICE, embodies an analytical ice accretion model 
that evaluates the thermodynamics of the freezing process that occurs when 
supercooled droplets impinge on a body. The atmospheric parameters of 
temperature, pressure, and velocity, and the meteorological parameters of 
liquid water content (LWC), droplet diameter, and relative humidity are 
specified and used to determine the shape of the ice accretion. The sur- 
face of the clean (uniced) geometry is defined by segments joining a set 
of discrete body coordinates (Fig. 2.1). The code consists of three major 
modules. They are l) the flow field calculation, 2) the particle trajectory 
and impingement calculation, and 3) the thermodynamic and ice accretion 
calculation. Each of these modules will be discussed in detail in following 
sections. 

LEWICE differs from other ice accretion prediction codes 1-4 because 
it applies a time-stepping procedure to “grow” the ice accretion. Initially, 
the flow field and droplet impingement characteristics are determined for 
the clean geometry. The ice growth rate on each segment defining the 
surface is then determined by applying the thermodynamic model. When 
a time increment is specified, this growth rate can be interpreted as an ice 
thickness and the body coordinates are adjusted to account for the accreted 
ice. This procedure is repeated, beginning with the calculation of the flow 
field about the iced geometry, then continued until the desired icing time 
has been reached. The application of this time-stepping procedure to the 
prediction of an ice accretion shape will be discussed in greater detail in 
following chapters. 

Ice accretion shapes for cylinders and several single-element airfoils have 
been calculated using this computer code. The calculated results have been 
compared to experimental ice accretion shapes obtained both in flight and 
in the Icing Research Tunnel at NASA Lewis Research Center. In general, 
the comparisons have been encouraging but inconsistent. Ice shapes for 
some conditions are well predicted, while for others there is little similarity 
between the predicted accretion and that obtained by experiment. Unfor- 
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tunately, the poorer predictions do not consistently occur in only one type 
of icing condition such as glaze or rime icing, on a specific type of airfoil, or 
in comparisons with results from a specific facility. There are many possible 
explanations for this behavior including inaccurate modeling of the accre- 
tion process, occasional errors in setting test conditions, and overextending 
the computer code into conditions where the assumptions, such as potential 
flow, do not apply. The known limitations will be addressed throughout 

this user’s manual. 

The manual begins with a discussion of the ice accretion process to 
identify the requirements of an ice accretion prediction methodology. The 
major components of the code and the mathematical models applied m 
each are then discussed. The input and output to the code are described, 
followed by the presentation of several sample cases, which illustrate various 
aspects of the code. The manual concludes with a summary of the current 
results and shortcomings of the method. 

Areas requiring additional development are identified and discussed in 
the Appendices (A through F). 




Figure 2.1: Geometry defined by segments joining body coordinates 
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Chapter 3 

DICUSSION OF THE ICE ACCRETION PROCESS 


The model of the ice accretion process applied in LEWICE is presented in 
this chapter, beginning with a discussion of some of the general character- 
istics of ice accretion shapes, and followed by a description of the physical 
model of the ice accretion process from which a mathematical model must 
be formulated. 

3.1 Ice Accretion Characteristics 

Before discussing the physical model of the ice accretion process, it is nec- 
essary to define some of the terms used in such a discussion. 

Ice may form on the forward facing surfaces of an aircraft flying through 
clouds composed of supercooled water droplets. The type and shape of ice 
that forms are functions of the atmospheric parameters of velocity, pres- 
sure, and temperature, and the meteorological parameters of liquid water 
content, droplet diameter, and icing time. 

Ice shapes are generally classified as glaze, mixed, and rime accretions. 
Rime ice is milky white and opaque, and will be denoted in the ice accretion 
profiles in this report by shading as shown in Figure 3.1a. Glaze ice is 
generally clear and is characterized by the presence of larger protuberances, 
commonly known as glaze horns, as shown in Figure 3.1b. A mixed ice 
accretion will have some of the characteristics of both glaze and rime ice 
accretions. As shown in Figure 3.1c, the center portion of a mixed ice 
accretion will have the characteristics of glaze ice accretion. This glaze 
center will be surrounded by rime ice accretions, commonly called rime 
feathers because of their thin, feather-like shape and delicate structure. 

The type of ice that will be formed is dependent on the atmospheric 
and meteorological conditions identified in the preceding paragraph. Pre- 
dicting the type and shape of the ice accretion that will be formed for a 
specified set of icing conditions is difficult because of the complex interac- 
tions between the atmospheric and meteorological parameters. Typically, 
rime ice is formed at lower temperatures, velocities, and LWC than glaze 
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ice. An example of the transition from glaze to rime ice as the total temper- 
ature is lowered is shown in Figure 3.2. At the warmest total temperature, 
Tt = —2.0 C, the accretion is composed exclusively of glaze ice. As the tem- 
perature decreases, areas of rime ice begin to form near the impingement 
limits. As the temperature decreases further, these rime portions increase 
in size until the accretion is composed solely of rime ice. The extent of icing 
and the locations at which ice forms on a surface are largely dictated by the 
size of the droplets impinging on the surface. For a given icing condition, 
and in the absence of ice shedding, the size of the accretion depends on 
the length of time ice is allowed to accrete. The general effects of temper- 
ature, droplet size, LWC, and angle of attack on ice shapes formed on a 
NACA 0012 airfoil in the NASA Lewis Icing Research Tunnel (IRT) are 
documented in Reference 6. 

3.2 Description of the Physical Model 

An understanding of the interactions between these parameters is required 
to predict the shape of an ice accretion that will be formed at a specified 
set of icing conditions. To develop this fundamental understanding, it is 
necessary to examine the physical model of the ice accretion process. 

A model of the ice accretion process was first presented by Tribus 7 and 
developed further by Messinger 8 . While many studies have been done to 
understand various aspects of the ice accretion process, the original physical 
model has been applied relatively unchanged. Recent close-up movies and 
photographs of the ice accretion process made at the NASA Lewis Research 
Center 9 have increased our understanding of the process and indicated that 
modifications to the physical model may be necessary. Conclusions drawn 
from the observations are included in the following discussion of the ice 
accretion process, and differences from the previous model are highlighted. 

The ice accretion process is characterized by the presence of supercooled 
droplets entrained in the flow about a body. These droplets follow trajecto- 
ries that will cause them to either be carried past or impinge upon a body. 
Upon impact with a clean surface, the droplets coalesce into larger surface 
drops under the effects of surface tension and flow along the surface as dic- 
tated by the airflow along the surface of the body. These surface drops will 
then either freeze on the surface or be shed from the surface because of the 


7 



aerodynamic forces on the drop. The ice accretions formed by this initial 
freezing form a rough surface which enhances the convective heat transfer 
and local collection efficiency of the surface, and therefore allows the ice 
accretion process to continue. 

The type of ice that will form for a given set of conditions is determined 
primarily by the rate at which the freezing process occurs. For example, if 
the conditions are such that the droplets freeze rapidly, there is essentially 
no initial coalescing and flowing of the droplets. Instead, they freeze on 
impact and form the characteristic rime ice accretions. These accretions 
are opaque and milky white in color because of the presence of air bubbles 
that are trapped in the structure during the rapid freezing process. As the 
rate of the freezing process decreases, the droplets begin to coalesce and 
flow on the surface. Upon freezing, these larger surface droplets form sur- 
face roughness elements which tend to enhance the convective heat transfer 
and local collection efficiency characteristics, which, in turn, enhance the 
continuing growth of the ice accretion in this region. These local areas 
of enhanced ice growth are, therefore, the beginnings of the characteristic 
horns found on mixed and glaze ice accretions. As the freezing rate de- 
creases further, the drops flow further along the surface of the body before 
freezing, thus moving the regions of enhanced ice growth away from the 
stagnation point. This, in turn, causes the horns of the accretion to move 
further apart and forms the familiar glaze ice accretions. As the rate of the 
freezing process decreases, less air is trapped within the ice structure and 
the ice gradually becomes clearer until it is essentially transparent, as in 
glaze ice. 

This description of the ice accretion process has identified four basic 
areas or processes that must be modeled in order to predict the shape of an 
ice accretion that will form on a body for a specified set of icing conditions. 
These areas are 1) the flow field about the body, 2) the droplet trajectory 
and impingement characteristics on the body, 3) the thermodynamics of 
the freezing process, and 4) the accumulation of ice on the surface of the 
body. To analytically predict an ice accretion shape, a mathematical model 
of each of these physical processes must be developed. Each of these areas 
will be discussed in subsequent sections to identify the methods used in the 
ice accretion prediction code, LEWICE. 


8 




Kim CONDITIONS 


_A_ 

JL 

JL 

CYLINDER diam, IN. 

1.0 

2.0 

1.0 

PSIA 

14.2 

12.2 

14.2 

V * 

-5.0 

21.1 

5 

v, fps 

200 

367 

400 

3m. u m 

20 

22 

15 

lwc, g/M 3 

0.6 

0.35 

0.53 

T, MIN 

8.0 

30.6 

6.8 


SHADED PORTIONS DENOTE RIME ICE — v 

\ 



a. Rime ice accretion, f = 1.0 



b. Glaze ice accretion, f = 0.15 


RIME FEATHERS- 
GLAZE 

center- 


's 



c. Mixed ice accretion, f = 0.53 

Figure 3.1: Examples of ice accretions formed at various atmospheric and 
meteorological conditions 
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(D) TOTAL TEMPERATURE, -18 °C 

(STAGNATION FREEZING FRACTION. 0.85). 


(E) TOTAL TEMPERATURE. -20 °C 
(STAGNATION FREEZING FRACTION, ) 


(F) TOTAL TEMPERATURE. -26 °C 

(STAGNATION FREEZING FRACTION. 0.9). 


Figure 3.2: Effect of temperature on the ice accretion shape. Thin ice 
samples removed from the airfoil and backlighted; Airspeed, 209km/hr, 
LWC, 1.3 g/m 3 ] DVM, 20m; Time, 8min; Airfoil, 0.53m chord 0012 airfoil 

at 4° angle. 


10 





Chapter 4 
METHODOLOGY 

The phenomena that make up the ice accretion process have been identified 
and now must be investigated in order to accurately predict the shape of 
an ice accretion that will form on a body. These include the evaluation of 

1) the flow field about the body, 

2) the droplet trajectory and impingement characteristics, 

3) the thermodynamics of the freezing process, and 

4) the accumulation of the ice on the surface. 

The computational methodology used to evaluate each of the above areas 
is discussed in the following sections. 

4.1 Calculation of the Flow Field 

The application of a flow field calculation method to an ice accretion pre- 
diction code requires that the method not only be able to calculate accurate 
flow fields about clean geometries, but also around the irregular, convoluted 
ice shapes that occur for many icing conditions. It is also desirable that the 
computational time and memory requirements of the method be as small 
as possible so that the use of the code is not limited to icing researchers 
with large computer facilities. 

The Douglas two-dimensional potential flow program developed by Hess 
and Smith, described in Reference 10, meets these requirements, and is used 
in the present study for calculating the flow field about the body. This pro- 
gram uses a distribution of sources, sinks, and/or vortices along the body 
surface to calculate the potential flow field. The body surface is represented 
by an arbitrary number of straight line segments. In calculating the flow 
field, contributions from all the sources, sinks, and/or vortices are summed. 
The accuracy of this method was tested by comparing its predicted veloc- 
ities and surface pressure coefficients with both analytical solutions and 
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experimental data 10 . Excellent agreement was found. Similar comparisons 
were performed at NASA Lewis and compared with pressure coefficient 
data found in Reference 11. These comparisons, shown in Figures 4.1a-j, 
show that the surface pressure coefficient is well predicted by the poten- 
tial flow method for angles of attack up to approximately 11.0 degrees and 
Mach numbers up to 0.50. These flow field limitations must be considered 
when predicting ice accretion shapes at high Mach numbers or high angles 
of attack. 

Only limited details of the methodology applied in the potential flow 
program are provided in this manual since the purpose of this study was to 
develop an analytical ice accretion prediction program. However, sufficient 
description of the input as well as subroutines are supplied to allow the user 
to become familiar with the primary aspects and limitations of the code. 

4.2 Calculation of the Droplet Impingement Characteristics 

The algorithm applied in the current ice accretion model to calculate the 
droplet trajectory and impingement characteristics was originally developed 
by Frost, Chang, Shieh, and Kimble of FWG Associates under contract to 
NASA Lewis 12-13 . 

This code was developed to calculate the trajectories and impingement 
characteristics of arbitrarily-shaped particles, and although water droplets 
are of primary concern when evaluating the ice accretion process, the gen- 
erality has been retained in LEWICE. While much of this code has been 
applied in its original form in the current study, there are significant changes 
in many areas and, therefore, the current version of the code is discussed 
in detail in this report. 

4.2.1 Definitions 

The primary droplet impingement characteristics that must be evaluated 
by an ice accretion prediction code are the regions of droplet impingement 
and the distribution of the mass of liquid on the surface of the body within 
this impingement region. 

The equations of particle motion, to be discussed in the following sec- 
tion, are used to calculate droplet trajectories from far upstream to impinge- 
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ment on the body and, if it occurs, the total and local collection efficiencies. 
The total collection efficiency is defined as the ratio of the actual mass of 
impinging water to the maximum value that would occur if the droplets 
followed straight-line trajectories. Figure 4.2 illustrates that this definition 
can be given in equation form as 


w 

where yo is the vertical distance between the droplet release points of the 
upper and lower surface tangent trajectories. The local collection efficiency, 
P, is also defined in Figure 4.2 and can be written in differential form as 


/? = 


dy o 
ds 


It is related to the total collection efficiency by the equation 


( 2 ) 

( 3 ) 


where s u and s t are the upper and lower surface impingement limits, respec- 
tively. The following sections will cover the methods applied to calculate 
the variables discussed above. 


4.2.2 Equations of Particle Motion 

The motion of a particle is analyzed as a point mass particle that is acted 
on by the potential flow field but which itself does not affect the flow. 
The forces acting on the particle are considered to be those of lift, drag, 
pitching moment, and gravity. Figure 4.3 shows the forces acting on the 
particle and the velocity vectors relative to the motion of the particle. 
The flight reference line (FRL) is not significant for a spherical particle; 
however, for arbitrarily shaped particles, i.e., a snow flake, the FRL must 
be defined relative to the lift, drag, and moment coefficient data available. 
The equations of motion of an arbitrarily shaped particle are derived from 
a force balance on a point mass, as shown in Figure 4.3, and are as follows: 

—* — ♦ t 

mi = —D cos 7 — L sin 7 + mg stn a 

my = —D sin 7 + L cos 7 — mg cos a (4) 
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where 


7 = tan 


-i 3/p Yjl 

x p — V x 


( 5 ) 


In Figure 4.3, note that the coordinate system used in LEWICE is fixed to 
the leading edge of the clean airfoil. 

For an airfoil at an angle of attack a, the coordinate system is at an 
angle to the gravitational coordinate system. Therefore, the effect of gravity 
must be accounted for in the equations for both lift and drag. 

The flow field velocity components in the x and y directions, i.e., V x 
and Vy, respectively, are obtained from the potential flow program. The 
aerodynamic drag and lift forces are defined as 


n * V * a 
D = c d —j- A p 

T V 2 A 

L = ci — — A p 


( 6 ) 


where A p is a characteristic area of the particle, p a is the density of air at 
the position of the particle, and V is the particle velocity relative to the 
flow field and defined as 

( 7 ) 


V = v/(ip - V.)’ + (if - V,)' 


For arbitrarily shaped particles, the pitch angle, 0 P , is required to eval- 
uate the angle of attack a p , using the following equation 

( 8 ) 


= Op ~ T 


This motion is governed by the following equation 


°=T 


( 9 ) 


where I tx is the moment of inertia of mass relative to the z axis. The 
moment of aerodynamic forces acting on the particle is 


M — c m A v d 


ip U ,p 


( 10 ) 
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where c m is the pitching moment coefficient which must also be specified 
by the user. 

The lift, drag, and pitching moment coefficients, cj, c d , and c m respec- 
tively, must be provided by the user for arbitrarily shaped particles. The 
coefficient data is input to the program through subroutine COEFF and 
should be functions of the particle angle of attack, as defined by Equation 
8, and the particle Reynolds number based on the particle diameter, given 
by the following equations: 

r., = v -± 

' v 

The diameter of the particle, d, and the kinematic viscosity of air, u, are 
assumed constant along the trajectory of the particle. 

Since water droplets are usually assumed to be rigid spheres in icing 
studies, the only forces considered to be acting on the particle are those of 
drag and gravity. The governing equations can therefore be simplified as 
follows: 

mi = —D cos 7 + mg sin a 

my = —D sin 7 — mg cos a (12) 

In this case, the drag force, D, is determined using a steady-state drag 
coefficient for a sphere which is a function of the droplet Reynolds number, 
Re p . Approximating droplets as rigid spheres is valid for drop radii less 
than 500.0 microns 14 . A valid drag law for spherical particles is built into 
the computer program in subroutine COEFF. 

For particles with diameters of less than 10 microns, the ratio of particle 
diameter to the mean distance between air molecules is small enough so that 
molecular slip phenomena result in drag forces lower than those calculated 
by the drag law used in LEWICE. The Cunningham correction factor, Cf, 
Reference 15, is therefore applied to correct the drag coefficient using the 
following equation: 

c d \slip = (13) 

G / 

The values of Cj are input by the user when necessary and are given in 
Table 4.1. As shown in Table 4.1, this effect is small for particles with 
diameters greater than 1.0 micron. Droplets this small would be included 
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Table 4.1: Cunningham Correction Factor for Standard Air 


d(u) 

C 

d(u) 

C 

0.001 

221.600 

0.1 

2.911 

0.002 

111.100 

0.2 

1.890 

0.003 

74.250 

0.3 

1.574 

0.004 

55.830 

0.4 

1.424 

0.005 

44.780 

0.5 

1.337 

0.006 

37.410 

0.6 

1.280 

0.007 

32.150 

0.7 

1.240 

0.008 

28.200 

0.8 

1.210 

0.009 

25.140 

0.9 

1.186 

0.010 

22.680 

1.0 

1.168 

0.020 

11.650 

2.0 

1.084 

0.030 

7.978 

3.0 

1.056 

0.040 

6.151 

4.0 

1.042 

0.050 

5.060 

5.0 

1.034 

0.060 

4.337 

6.0 

1.028 

0.070 

3.823 

7.0 

1.024 

0.080 

3.441 

8.0 

1.021 

0.090 

3.145 

9.0 

1.019 



10.0 

1.017 


in an ice accretion prediction only in exceptional circumstances, such as 
testing the limiting capabilities of an ice accretion code. 


4.2.3 Method of Integration 

The equations governing the motion of arbitrarily-shaped particles are as 
follows: 


x = 


dx 

dt 



e 

* = d * 


dx 

~dt 


D L . . 

cos 7 sin 7 + g sin a 

m m 


(14a — c) 
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dy D . L 

— = sin 7 H cos 7 — g cos a 

at m m 

dti_ _ M_ 
dt I zz 

(14 d-f) 


For spherical water droplets, Equations 14c and 14f are not applicable 
and Equations 14d and 14e are simplified to result in the following four 
equations: 


x — 


dx 

dt 



dx D , 

— cos 7 + g sin a 

dt m 

dy D . 

— = sin 7 — g sin a 

dt m 


(15) 


These equations are integrated using the method of Gear developed for stiff 
equations 16-17 . The details of the subroutines that make up the integration 
method, i.e., DIFSUB, DECOMP, SOLVE, and PEDERV, can be found in 
Reference 16 and in COMMENT statements in the computer code. The 
integration routine also requires that the equations to be integrated be 
located in a subroutine named DIFFUN. This subroutine currently contains 
Equations 14 - 15. 


4.2.4 Determination of Droplet Impingement 

The calculation of the droplet trajectories is continued until the droplets 
impinge upon the body or move out of range. This section describes the 
procedure used to determine whether or not a droplet impacts the body 
and, if so, the location of impingement. These calculations are controlled 
by subroutine MODE. 
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As previously discussed, the geometry is defined by segments joining a 
discrete set of body coordinates as shown in Figure 2.1. A droplet is con- 
sidered to impact the body when its trajectory intersects one of these body 
segments. The current model does not take account of grazing collisions 
or droplets that may impact the body so that they are re-introduced into 
the flow by bouncing, splashing, etc. The impact algorithm, found in sub- 
routine INTRST, sequentially sums the angles between lines drawn from 
the particle position to adjacent points describing the closed curve of the 
geometry, as shown in Figure 4.4. The summation starts with the angle 
between lines drawn to the first and second points, continuing with the 
angle between the lines to the second and third points, and so on, all the 
way to the angle between lines drawn to the next to last and last points. 
If the particle is outside the body, the sum of these angles will always be 
zero. If the particle has crossed one of the body segments and lies inside 
the body, the sum of the angles is 2n. If the particle lies directly on one of 
the body segments, the summed angle will equal 7 r. 

A particle trajectory is calculated until the summed angles total n or 2n, 
which indicates that the particle has impinged upon the body, or until the 
particle passes outside of the pre-specified boundaries. The impact point 
results from the intersection of the particle trajectory and the line con- 
necting the adjacent ice shape points through which the particle trajectory 
passed (Figure 4.5). The particular line segment through which the particle 
passed is determined by first calculating the intersection of the line joining 
the present and previous particle positions and the line formed by each of 
the adjacent points that describe the body geometry, as shown in Figure 
4.5. If the particle passed through a particular segment, the distance from 
the intersection (impingement) point to the endpoints of the segments will 
be less than the length of either the trajectory or body segments. Once the 
body segment through which the particle passed and the intersection (im- 
pingement) point have been determined, the surface distance, s, from the 
stagnation point to the impingement point is determined by interpolation. 

4.2.5 Calculation of the Local Collection Efficiency 

The particle trajectories and impingement points, calculated as previously 
described, are used to establish the relations between the particle’s initial 
position (x 0 ,yo) and the position where it impinges on the body surface, 
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specified by the surface distance, s, which is the length along the body 
surface measured from the stagnation point. The value of s is defined as 
negative on the lower surface and positive on the upper surface. 

The local collection efficiency is calculated by first calculating droplet 
trajectories and producing a plot of particle release point, yo, vs. surface 
impact distance, s, as shown in Figure 4.6. As was indicated by Equation 
2, the local collection efficiency is a function of the surface distance and 
can be determined by differentiating the curve shown in Figure 4.6 with 
respect to s. 

The derivative at the center of each body segment is calculated by first 
determining the four yo vs. s points whose s values are closest to the s 
value of the body segment at which the local collection efficiency is desired, 
as shown in Figure 4.6. These four points are then fit with a quadratic 
polynomial using the method of least squares. The local collection efficiency 
at the desired s location is determined by differentiating the polynomial. 
The local collection efficiency is calculated in subroutine EFFICY, while 
the curve fitting/differentiation procedure is found in subroutine TERP. 

4. 2. 5.1 Local Collection Efficiency Calculation for Multidispersed 
Particle Distributions 

The previous section described how the local collection efficiency was calcu- 
lated for a single droplet diameter. In icing applications, the mass median 
droplet diameter of the droplet size distribution is used to characterize the 
size of the droplets. A feature of the particle trajectory portion of the tra- 
jectory program is that it allows the user to analyze the local collection 
efficiency for a multidispersed particle distribution. 

To perform this calculation, the user must input the droplet diameter 
and the associated mass fraction and Cunningham correction factor for 
each specified droplet size. A maximum of 10 droplet sizes can be used to 
characterize a droplet distribution. For example, the required input for a 
Langmuir D distribution with a mass median of 20.0 microns is shown in 
Table 4.2. 

The solution procedure is begun by calculating the local collection effi- 
ciency distribution for each droplet size characterizing the distribution. The 
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Table 4.2: Langmuir D droplet size distribution with a mass median of 20 
microns 


Percentage 

Ratio of 

LWC 

Diameters 

0.05 

0.31 

0.10 

0.52 

0.20 

0.71 

0.30 

1.00 

0.20 

1.37 

0.10 

1.74 

0.05 

2.22 


Droplet 
Diameter (M m ) 
6.2 

10.4 
14.2 
20.0 

27.4 
34.8 

44.4 


Cunningham 
Correction Factor 
1.0272 

1.00 

1.00 

1.00 

1.00 

1.00 

1.00 


local collection efficiency for the distribution is determined by summing the 
contributions of each of the droplet sizes using the following equation: 

P (s) = J2 n< ft ( s ) ( 16 ) 

•=i 

where n, is the mass fraction of liquid water associated with droplet di- 
ameter i and N is the number of droplet sizes used to characterize the 
distribution. The local collection efficiency for a droplet size distribution is 
also calculated in subroutine EFFICY. 

4.2.6 Computational Procedure 

The previous sections described various aspects of the calculation of the 
particle impingement characteristics. The purpose of this section is to 
describe how these calculations work together to yield the desired local 
collection efficiency information. 

After calculating the flow field about the body, the program enters the 
particle trajectory main subroutine, TRAJ. First, the initial particle loca- 
tion *o, yo and velocity V Xp ,V Vp must be determined, either by the computer 
program or from information input by the user. A particle should be re- 
leased at a location upstream of the airfoil where the flowfield is essentially 


20 


the same as the free stream conditions. The program will select an ini- 
tial upstream x-coordinate, i 0 » by searching for a position where the local 
velocity V L and the freestream velocity Voo satisfy the following inequality: 


1 - 


V L 


'00 


Vo In 

VoU 


< € 


(17) 


where e (VEPS in the computer program) is specified by the user. Equation 
(17) is tested over a specified range, y 0 | m , n < y a < J/olma*. where y Q \min 
and t/o|mazt illustrated in Figure 4.7, are also input by the user. 

With the initial upstream x-coordinate known, the next step is to locate 
two trajectories, one that passes above the body and one that passes below 
the body. The vertical distances from which these particles are released, 
yolmaz and t/o | m i n , respectively, are specified by the user. The calculation of 
these particle trajectories is controlled by subroutine RANGE. Using these 
upper and lower trajectories as boundaries, the upper and lower impinge- 
ment limits, you and yoi, are determined in subroutine IMPLIM (Figure 
4.7b). This subroutine uses a Newton iteration scheme to determine the 
release points for particles with trajectories that impinge upon the body 
tangent to the surface. For example, when searching for the upper impinge- 
ment limit, the trajectory of a particle released from j/o|i = , s 

computed. If the particle passes under or hits the body, the next trajectory 
is computed from y 0 |2 = troh + r° l wtil . if it passes over the body, the next 
trajectory is calculated from y 0 |j = (y°h+|[ 0 l n» id. Successive halving of the 
range y 0 | m m to yo|moz continues until the upper impingement limit is found. 
Convergence of this iterative procedure is assumed when the difference be- 
tween the y 0 values of two trajectories, i.e., one that hit the body and one 
that missed, is less than a small value specified by the user (YOLIM). This 
procedure is then repeated for the lower surface to determine the lower 
surface impingement limit. Any particle released between the upper and 
lower impingement limits will strike the body and any particle released 
from outside this range will miss the body. 

With the limiting release points known, the program enters subroutine 
COLLEC where the range of vertical position, yo u to y 0 j, is divided into a 
number (NPL) of equally-spaced increments prescribed by the user. The 
trajectory of particles leaving each of these vertical positions is calculated, 
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and the impingement position of the particle on the body surface is deter- 
mined using the method previously described. The y 0 vs. surface distance, 
s, information obtained in this calculation is then differentiated in subrou- 
tine EFFICY to determine the local collection efficiency at the midpoint of 
each of the body segments. 

4.3 Calculation of the Thermodynamic Characteristics 

The thermodynamic analysis of an icing surface was first developed by 
Tribus 7 from the physical model of the ice accretion process previously 
discussed. This model was used to calculate the heating requirements for 
icing protection and proposed LWC measurement systems. Messinger 8 de- 
veloped the thermodynamic model further to include an analysis of the 
temperature of an unheated surface in icing conditions for three surface 
temperature regimes, i.e., less than 273.15 K, equal to 273.15 K, and above 
273.15 K, and the concept of the freezing fraction, f, to be discussed later. 
These early formulations have been used in various icing applications. 

As discussed in Section 3.2, microscopic movies of the ice accretion pro- 
cess made at NASA Lewis Research Center® indicate that the process may 
be more accurately modeled by modifying the equations used in past icing 
studies. The observations reveal that, after the initial flow of the coa- 
lesced droplets on the surface, the liquid does not flow but is caught and 
frozen in the grooves between the individual surface roughness elements. 
Incorporating this observation into a mathematical model would proba- 
bly require modeling the individual roughness elements and the freezing of 
pools of water surrounded on all sides by ice. A microscopic and possibly 
three-dimensional analysis of the icing surface would be required to math- 
ematically apply this model to an ice accretion prediction method. The 
mathematical model used in this and previous studies is more macroscopic 
in nature because the roughness elements do not directly effect the freezing 
process except to enhance the convective heat transfer coefficient. 

The equations that model the thermodynamics of the freezing process 
on a body undergoing icing are formulated by performing a First Law of 
Thermodynamic mass and energy balance on a control volume located on 
the surface. The control volume to be analyzed is located on the surface of 
the body and extends from outside the boundary layer to the surface of the 
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body, as shown in Figure 4.8a. The lower boundary of the control volume is 
initially on the surface of the clean geometry and moves outward with the 
surface as the ice accretes. Therefore, the control volume is always situated 
on either the clean or ice surface. Computationally, a control volume is 
placed over each segment defining the body geometry, as shown in Figure 
4.8b. The equations resulting from the mass and energy balance can be 
expressed as follows: 

Mass Balance: 


= (1.0 - /)(m e + m r .J - m e (18) 

Energy Balance: 

m c [w(r, -273.15) + ^f] + 

rh fin [c pW)4U r(,-x)(r 4ur (, +1 ) - 273.15)] + q k As = 

*^e[c pW( * ur [(T, ur — 273.15) + Z»y] + 

[(1 - f)(th e + m r ,.J - m e ]c pWitur (T tur - 273.15)+ 
fm e - m r; J [cp t -, Jur (r tur - 273.15) - L/] + 

h e \T. ur -T L - r -^£\As (19) 

*Cpa J 

The complete derivation of Equations (18) and (19) is included in Appendix 

A. 


4.3.1 Definitions 

Before discussing the method used to solve Equations (18) and (19), it is 
necessary to discuss several of the terms that are important to the solution 
procedure. 

As discussed in Section 3.1, the atmospheric and meterological param- 
eters determine the type of ice that will form for a given icing condition. It 
has been found by various authors that the concept of a freezing fraction 
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can be used to determine the type of ice that will form. The freezing frac- 
tion, /, was defined by Messinger as the fraction of impinging liquid that 
freezes within the region of impingement. In this application, f is defined 
as the fraction of the total liquid entering the control volume that freezes 
within the control volume. It is given by the equation 


/ = 


m, 

the + th rin 


( 20 ) 


For colder icing conditions, the droplets tend to freeze immediately on 
impact, resulting in the formation of rime ice. Since all the water entering 
the control volume freezes within the control volume, the freezing fraction 
equals 1.0. Freezing fractions close to 0.0 characterize glaze or clear ice. 
Freezing fractions between approximately 0.3 and 1.0 will normally indicate 
that the ice has some combination of glaze and rime characteristics. As 
shown in Figure 3.1, ice accretions are often composed of glaze, rime, and 
intermediate regions. The local value of the freezing fraction therefore 
varies along the surface, and can be calculated using the mass and energy 
balances given by Equations (18) and (19). 


4.3.2 Solution of the Energy Equation 

The evaluation of Equation (19) is begun at the stagnation point because 
there will be no runback into the control volumes located on each side of 
the stagnation point, as shown in Figure 4.8b. Therefore, 

m rin \ttag — rh r| .J,j og _i = 0.0 (21) 

It is first assumed that the equilibrium surface temperature, T >ur , equals 
273.15 K. The terms of Equation (19) are then evaluated at this tempera- 
ture, and the resulting expression is solved to determine the freezing frac- 
tion, /. This calculation is performed in subroutine COMPF. The value 
of / will be either 1) less than 0.0, 2) between 0.0 and 1.0, inclusive, or 3) 
greater than 1.0. 

For 0.0 </ < 1.0, T lxir = 273.15#, and the initial assumption was 
correct. A value of / < 0.0 indicates that the surface temperature is greater 
than 273.15 K. Therefore, the solution is obtained by setting / = 0.0 and 
solving for T, ur in subroutine COMPT. Note that an iterative procedure 
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is required since many of the terms of Equation (19) are functions of T tur . 
Similarly, / > 1.0 indicates that T tur is less than 273. 15K, and / should be 
set equal to 1.0. Again, an iterative procedure must be applied to determine 

T, ur . 

When the thermodynamic characteristics of the control volume are 
known, the mass balance given by Equation (18) is used to determine the 
mass flow rate of runback water out of the control volume. Any water flow 
out of the control volume will be away from the stagnation point and into 
the next control volume. 

The above procedure is then repeated for the adjacent downstream con- 
trol volume and continued along the upper surface of the body. The entire 
procedure is then repeated again, starting at the stagnation point and pro- 
ceding along the lower surface of the body. 


4.4 Calculation of the Iced Geometry 

When the freezing fraction has been determined for each segment (con- 
trol volume) on the body, Equation (20) is used to calculate the local ice 
accumulation rate, rewritten below as 

mi = f{m e + m r ,.J (22) 

This ice growth rate must be interpreted as an ice thickness to form an 
ice accretion on the surface of the geometry. The thickness of the ice layer 
grown on a particular segment is given by the equation 

4 = — — — — — — (23) 

Pi 

where />,• is the density of the ice, As is the length of the segment, and At 
is a time increment specified by the user. 

The density of the accreted ice is determined using the empirical expres- 
sion developed by Macklin (Reference 18). This correlation was developed 
from predominantly rime ice accretions at low temperatures and velocities 
and is as follows: 
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(24) 


Pi 



y 76 

2 T IU r ’ 


In this expression, d m is the mass median droplet diameter in microns, V d 
is the droplet impact velocity in m/sec, and T, ur is the surface temperature 
in °C. In the calculation, the freestream velocity V*, is used for V d . The ice 
density has the units of kg/m 3 . Equation (24) is used to determine the ice 
density when the atmospheric and meteorological parameters are such that 

pm m < -d m Vj < pm m 

sec °C ~ 2 T, ur ~ sec°C 

and T. ur < -5 °C. When these conditions are not satified, it is assumed 
that the ice has a density of 917 kg/m 3 . In general, the inequality will be 
satisfied under conditions of small droplets, low velocities, and low surface 
temperatures. For example, with a droplet diameter of 12 microns and 
a velocity of 60 m/s, the surface temperature must be less than -20 °C, a 
rather extreme rime condition. The density of the ice accretion is evaluated 
for each surface segment to allow mixed ice accretions to form (mixed ice 
accretions contain both rime and glaze ice). The calculation of the ice 
density is performed in function RHOICE. 

The new ice surface is formed by first adding the ice thickness, d,, per- 
pendicular to each segment, as shown in Figure 4.9. The adjacent endpoints 
of each of these new segments are then averaged to obtain the coordinates 
describing the new ice surface. The new ice surface is calculated in subrou- 
tine NWFOIL using this procedure. 

When the new surface is formed, the length of a segment increases. 
The segments are allowed to grow but, at some point, must be split to 
maintain adequate definition of the surface. The segment length is allowed 
to increase until it is SEGTOL times the length of the original segment. 
The value of SEGTOL is input by the user. When a segment has grown 
to SEGTOL times the initial segment length, it is divided in half to form 
two new segments, and a new point is added to the set of coordinates. 
These two new segments will be allowed to grow to SEGTOL times their 
current length before being divided. Note that a value of SEGTOL < 2.0 
will result in progressively shorter segments. Values greater than 2.0 will 
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result in progressively longer segments. The segment lengths are checked 
and, if necessary, divided in subroutine NWPTS. 

As the ice accretion grows, it is also possible for two lobes of the accre- 
tion to grow together, causing some of the points to lie in the interior of 
the body, as shown in Figure 4.10. These points must be removed in order 
to continue the calculations. 

The point removal procedure is begun by applying the same procedure 
used to determine if two segments intersect, as shown in Figure 4. In this 
case, each body segment is checked with every other segment, excluding the 
two adjacent segments. As shown in Figure 10, if an intersection is found, 
all segments between the two intersecting segments are removed, and the 
set of coordinates is revised to reflect these changes. This procedure is 
performed in subroutine SEGSEC. 

4.5 General Computational Procedure 

The previous sections discussed each of the individual phenomena of the 
ice accretion process that are evaluated in LEWICE. The purpose of this 
section is to describe how these individual calculations are implemented to 
form a complete ice accretion. 

As discussed in the Introduction, LEWICE applies a time-stepping pro- 
cedure to grow an ice accretion. The flow field and droplet impingement 
characteristics are initially determined for the clean geometry. The ice 
growth rate on each segment defining the surface is determined by applying 
the thermodynamic model. The new surface is then formed by specifying 
an icing time and applying the procedure described in the previous section 
to account for the accreted ice on the clean surface. After calculating this 
initial ice layer, two options are available. 

The most desirable option is to repeat the entire procedure, beginning 
with the calculation of the flow field about the iced geometry, to obtain 
revised local collection efficiency and thermodynamic data. Unfortunately, 
since the majority of the computational time is spent calculating droplet 
trajectories, this option also increases the computational time required to 
accrete a layer of ice. Therefore, a second option was made available. 
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If the amount of ice accreted during the time step is small or no new 
protuberances, such as glaze horns, were formed, it is possible to accrete 
another layer of ice using the same local collection efficiency curve calcu- 
lated from the previous time step. This option, of course, does not produce 
results that are as accurate as those in the first option, especially for glaze 
ice accretions. The advantage is that the computational time required is 
significantly reduced. 

Each of these options will require that another time increment be spec- 
ified. The above procedure is repeated by specifying discrete time incre- 
ments until the desired icing time is reached. Guidelines for choosing an 
appropriate time increment are also given in Section 5.3.1. 

Ice accretions can have many geometrical shapes ranging from the 
smooth, aerodynamically-shaped rime accretions to rough glaze accretions 
with deep center grooves. LEWICE is therefore required to calculate suf- 
ficiently accurate flow fields and particle trajectories about what can be 
very irregular geometries where viscous effects such as boundary layer sep- 
aration and reattachment are important. This can, at times, exceed the 
capabilities of the potential flow code and produce non-physical results. 
Much work has been done to identify and correct inaccuracies in the flow 
field calculations. The techniques that have been implemented in the code 
to overcome these flow field inadequacies will be discussed in many of the 
following sections. 
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Figure 4.1: Continued. 





s u = Upper — Surface Impingement Limit 
si = Lower — Surface Impingement Limit 
H = Forward Projection of the Airfoil Height 

Total Collection Efficiency 


Local Collection Efficiency 



ds 


Figure 4.2: Definition of total and local collection efficiency 



Figure 4.3: Forces acting on an arbitrarily-shaped particle 
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N = number of points describing the geometry (point 1 is the same as 
poing N). 

If z p , t/p lies outside the body, 0 T = 0. If x p , y p lies on the body, 0 T = n. If 
Zp, y p lies inside the body, Or = 27T 

Figure 4.4: Illustration of the method to determine particle impact 


r- PARTICLE POSITION 
I AT PREVIOUS 




I. Particle does not pass through II. Particle passes through the 

the segment being evaluated. segment being evaluated. 

The particle passed through the segment being evaluated if all of the fol- 
lowing criteria are satisfied: 


IPi < PoPi 


JN i+l < N< N i+l 

7F 0 <PoP 1 

W < Ni N i+l 

Figure 4.5: Illustration of the method to determine the location of particle 
impingement 
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Figure 4.6: Particle release position, yo, vs. surface impace distances, s 
points used to determine the polynomial 
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a. Criteria to determine the particle release location. 



b. Particle release locations for the upper and lower surface tangent trajectories 
Figure 4.7: Definition of terms used to determine the Local Collection Efficiency. 
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a. Single control volume on the icing surface. 



b. Thermodynamic control volumes over each segment defining the body geometry 

Figure 4.8: Identification of the control volume used to formulate the ther- 
modynamic equations 
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Figure 4.9: Illustration of the method to calculate the iced geometry 
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a. Geometry with intersecting segments (5)and(l) 



b. Revised geometry with the intersecting segments removed. (The original 
segment numbers from a. axe in parenthesis) 

Figure 4.10: Illustration of the method to remove points from the body 
geometry 



Chapter 5 

LEWICE INPUTS 


This chapter presents the input format for a LEWICE input card deck. 
Section 5.1 contains the documentation that is needed to set up, execute, 
and make changes to the data deck. When more explanation is required, 
the user should consult Section 5.2 which includes additional information 
concerning input definitions and program options. Examples of the input 
files are shown in Section 7.0. The interactive input requested by the pro- 
gram is discussed in Section 5.3. 


5.1 Input Format 

The purpose of this section is to provide a quick reference to the input 
parameters used in LEWICE. Additional information can be found in Sec- 
tion 5.2 and in various sections throughout the text. Where applicable, 
the sections that contain additional information about the parameters are 
identified. 

5.1.1 Potential Flow Input 

CARD 01 Run Identification Card (8A4) 

IDR 

CARD 02 Potential Flow Code Control Parameters (NAMELIST S24Y) 

Variable Description 

Lift Control Flag 
= 0 This is not a lifting body 
= 1 This is a lifting body 


ILIFT 



IPARA 

IFIRST 

ISECND 

IPVOR 

INCLT 

CLT 

ICHORD 

CCL 


Element Geometry Flag 
= 0 Linear Elements 
= 1 Parabolic Elements 
First-order Terms Flag 
= 0 No first-order terms 
= 1 First derivative term 
= 2 Curvature term 
= 3 Both first-order terms 
Second-order terms flag 
= 0 No second- order terms 
= 1 Second derivative term 
= 2 Curvature squared term 
= 3 Both second- order terms 
Vorticity Distribution Flag 
= 0 Use constant vorticity between 
body elements 
= 1 Use variable vorticity 
distribution 
cj, a Flag 

= 0 Angle of attack, a, is input 
= 1 Total lift coefficient, cj, is input 
Value of angle of attack (degrees) 
or lift coefficient depending on the 
value of INCLT 
Reference Length Flag 
= 0 The reference length used in 
calculating c< is to be set = 1.0 
= 1 The reference length used in 

calculating c j will be input as CCL 
The value for the reference length 
(chord) used in calculating c j 
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IND 


ISOL 


IPRINT 


IFILL 


Individual Solution Flag 
S24Y is capable of calculating the 
potential flow about up to 6 bodies 
and then superimpose the results of 
each. The possible values of IND are 
as follows: 

= 0 Edge velocities are not calculated 
for each body 

= 1 Edge velocities are calculated 
for each body 

In LEWICE, only one body is input and the 
edge velocities are always required. 

Therefore, IND = 1 

Matrix Solution Method Control Flag 

= 0 Use routine SOLVIT for the matrix 

solution (used when a very large number 
of geometry points have been input) 

= 1 Use routine QUASI for the matrix 
solution 

= 2 Use routine MIS1 for the matrix 

solution. Maximum number of geometry 
points is 101. If the number of 
points is greater than 101, the program 
will automatically use SOLVIT. 
Print/Punch Flag 
= 0 Normal output 
= 2 Print the individual matrices 
= 7 Punch the output on cards 
IPRINT should be set equal to 0 to 
reduce the amount of printed output 
Parabolic Integration Flag 
S24Y calculates the forces and moments 
acting on the body using both trapezoidal 
and parabolic integration of the calculated 
pressure coefficient. The results of the 
trapezoidal calculations are always output. 
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The value of IFILL determines whether the 
parabolic results are printed. 

= 0 Results of the parabolic integration are 
not printed 

= 1 Results of the parabolic integrations are 
printed 

ICOMB Combination Solution Flag 

= 0 No combination solution calculated 
= 1 Combination solution calculated 


CARD 03 x-Coordinates (6F12.7,2X,I1,1X,I1,1X,I1) 


Column 

Variable 

Description 

01-12 

X(1) 

x-coordinate of the geometry. Up to 

13-24 

X(2) 

six coordinates may be input on each 

25-36 

X(3) 

card depending on how the INO flag i 

37-48 

X(4) 

set. 

49-60 

X(5) 


61-72 

X(6) 


75 

INO 

Number of data points per card. If 
there are 6 values per card, INO may 
be left blank. 

77 

ISTAT 

Last Card Flag 
=0 This is not the last 
x-coordinate card. More cards 
will follow. 

=1 This is the last x-coordinate card. 

79 

ITYPE 

x-Coordinate Flag 
= 3 
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CARD 04 y-Coordinates (6F12.7,2X,I1,1X,I1,1X,I1) 


Column 

Variable 

Description 

01-12 

Y(l) 

y-coordinate of the geometry. Up 

13-24 

Y(2) 

to six points may be input on each 

25-36 

Y(3) 

card depending upon how the INO flag 

37-48 

Y(4) 

is set. 

49-60 

Y(5) 


61-72 

Y(6) 


75 

INO 

Number of data points per card. If 
there are six values per card, INO 
may be left blank. 

77 

ISTAT 

Last Card Flag 

= 0 This is not the last y-coordinate 
card. More cards will follow. 

= 1 This is the last y-coordinate 
card. 

79 

ITYPE 

y-Coordinate Flag 
= 4 


5.1.2 Trajectory Code Input 

CARD 05 Tajectory Input I (NAMELIST TRAJ1) 


Variable 

Description 

GEPS 

Convergence criterion for the integration 


method of Gear 

VEPS 

Accuracy criterion for the case when 


LXOR = l (Section 4.2.6) 
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DSHIFT 

LCMB 

LCMP 

LEQM 

LSYM 

LYOR 

LXOR 


x-distance to shift coordinates after 
the potential flow calculation to avoid 
discretization errors 
Combination Correction Flag 
= 0 Calculates the combination solution 
using the method of S24Y 
= 1 Calculates the combination solution 
using the method of COMBIN-2D 
Compressiblity Correction Flag 
= 0 No correction for compressibility 
= 1 Correct velocity values to account 
for compressibility 
Particle Initial Condition Flag 
= 0 The initial x- and y-components of 
the particle velocity will be input 
= 1 The initial particle velocity is 
equal to the flow at the initial 
particle location (in equilibrium 
with the flow) 

Symmetric Flow Field Flag 
= 0 Unsymmetric flow field (general case) 

= 1 Symmetric flow field (only half plane 
is computed) 

y-Coordinate Particle Release Flag 
= 0 Particle is released from the position 
specified by YORC 

= 1 Program determines the vertical particle 

release position using YOMAX and YOMIN as 
the initial guesses 
x-Coordinate Particle Release Flag 
= 0 Particle is released from the position 
specified by XORC 

= 1 Particle release position is determined 
using the criteria |l — < V EPS 
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NEQ Number of equations to be solved to 

determine the particle trajectories 
= 4 Spherical, non-lifting particle 
= 6 Lifting, rotating particle 
NPL Number of particle trajectories to be 

computed to define the yo vs. s curve. 

If NPL is set equal to 1, a single 
trajectory is to be calculated. 

NSEAR Maximum number of trajectories allowed 
to be calculated in the search for the 
upper and lower impingement limits 
NSI Number of droplet sizes used to 

characterize the cloud droplet 
distribution (maximum of 10) 

TIMSTP Initial value of the time step used in 
the integration of the particle 
trajectory equation (Gear’s integration 
method) 

CARD 06 Trajectory Input II (NAMELIST TRAJ2) 


Variable 

Description 

CHORD 

Airfoil chord (m) 

G 

Acceleration of gravity (m/s 2 ) 

PIT 

Initial angle of the particle flight 
reference line (Figure 4.3) (degrees) 

PRATK 

Initial value of the particle angle 
of attack (Figure 4.3) (degrees) 

XORC 

x-coordinate position of particle 
release (xo/chord). XORC need not 
be input if LXOR = 1. 

YORC 

y-coordinate position of particle 
release (yo /chord). YORC need not be 
input if LYOR = 1. 
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XSTOP 


Maximum downstream distance, normalized 
by the chord of the airfoil, for which 
particle trajectories are calculated 
YOLIM Accuracy criterion for computing the 
surface impingement limits 
(Section 4.2.6) 

YOMAX Initial guess for the y-coordinate of 
the upper surface tangent trajectory 
release point (yo/chord) 

(Section 4.2.6) 

YOMIN Initial guess for the y-coordinate of 
the lower surface tangent trajectory 
release point (y 0 /chord) 

(Section 4.2.6) 

VXPIN, x- and y-components of the particle 
VYPIN release velocity. They are input only 
when LEQM = 0. 

CARD 07 Droplet Distribution Characterization Card 
(NAMELIST DIST) 

Variable Description 

DPD Droplet sizes in the distribution 

(maximum of 10) (microns) 

FLWC Fraction of LWC for each droplet size 
specified in the distribution 
CFP Cunningham correction factor 

(Section 4.2.2, Table 4.2) 


5.1.3 Ice Accretion Input 

CARD 08 Ice Accretion Data (NAMELIST ICE) 
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Variable 

Description 

VINF 

Free-stream velocity (m/s) 

TAMB 

Static temperature (K) 

PAMB 

Static pressure (Pa) 

LWC 

Liquid water content (g/m 3 ) 

DPMM 

Mass median droplet diameter of the 
specified droplet distribution (microns) 

RH 

Percent relative humidity 

XKINIT 

Initial value of the equivalent sand- 
grain roughness of the icing surface 
(Chapter 4.3, Appendix F) 

SEGTOL 

Maximum amount any segment may grow 
before it is divided in two 
(Chapter 4.4) 


5.2 User’s Guide to Input Format 
5.2.1 Potential Flow Input 

As discussed in Section 4.1, the flow field used in LEWICE is calculated 
using the Douglas potential flow method (S24Y). This code was developed 
to calculate the flow field about a wide variety of geometries. Using this 
code in an ice accretion prediction method decreases the generality required 
by the potential flow code, and many of the input parameters and program 
options are not necessary for this application. The input to the potential 
flow code required by LEWICE users was simplified to include only the 
applicable parameters. However, the potential flow computer code was 
not modified and remains in its original form in LEWICE. Therefore, the 
original input format to the potential flow code is still used and all input 
parameters are required. The input parameters not found in NAMELIST 
S24Y are set to default values in subroutine SETUP. In this subroutine, 
the input file to the potential flow code is set up and written to unit 45. A 
description of the complete input to the potential flow code written to unit 
45 is given in Appendix C. Further details concerning the input parameters 
can be found in Reference C-l and C-2. 
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5.2.2 Trajectory Code Input 

This section provides the user with additional information about many of 
the input parameters to the particle trajectory code. Suggested values of 
the parameters are also given. 

5. 2. 2.1 Trajectory Input 1 (NAMELIST TRAJ1) 

GEPS 

This variable is the error test constant used in the integration method 
of Gear. Single step error estimates made in the integration algorithm must 
be less than GEPS in the Euclidean norm. The step size and/or order is 
adjusted so that this criteria is met. See References 16 and 17 for additional 
information on the parameter and the integration scheme in general. 

Parametric studies were made to evaluate the effect of changing the 
value of GEPS from 0.001 to 0.00001. Increasing the value of GEPS was 
found to reduce the cpu time required to calculate each trajectory by allow- 
ing larger integration step sizes to be used. Figure 5.1 shows how increasing 
GEPS decreases the computational time per trajectory. 

Unfortunately, larger values of GEPS also allow larger computational 
errors which are reflected in the calculated particle impingement locations. 
Computationally, GEPS should approach 0.0, but the required compu- 
tational time would be excessive, as shown in Figure 5.1. A value of 
GEPS=0.00005 was therefore used for all calculations. 

VEPS 

This parameter is used when the program is to locate the proper x 
location at which to release the particles. The particles will be released 
from an x location where, between Y0MIN and Y0MAX, 

|l - ^1 < VEPS 

*oo 

Increasing the value of VEPS will allow particles to be released closer to 
the body and thereby decrease the number of integration timesteps required 
for the particle to strike the body. The computational time will therefore 
be decreased. 
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No studies have been made to show the effect of releasing the particles 
closer to the body. The particles must be released in essentially free stream 
conditions for physically accurate trajectories to be calculated. For this 
reason, a VEPS value of 1.0 x 10~ s has been used for all calculations. 

DSHIFT 

The potential flow computer program has a relatively large discretiza- 
tion error very close to the body. Figure 5.2, taken from Reference 12, 
shows the longitudinal and vertical velocities around the leading edge of a 
Joukowski airfoil. The velocity is computed for different constant values 
of separation distance (DSHIFT) from the body, as illustrated in the in- 
sert. Note the large oscillations in the flow field velocity near the surface. 
These oscillations can cause erratic particle trajectories close to the body, 
especially for small particles that are effected by small flow field perturba- 
tions, and can cause fatal program errors to occur, thereby terminating the 
program. 

To overcome the effect of the discretization error near the body, an ar- 
tificial impingement surface is generated by the computer program. This 
surface is defined by displacing each point of the body by a small incre- 
ment DSHIFT in the upstream x direction. This displacement essentially 
increases the size of the body to include the region where discretization 
errors are present. DSHIFT values of approximately .2 to .6 percent of the 
chord length are commonly used. If irregularties in the impingement curves 
or droplet trajectory calculations persist for a specific case, the DSHIFT 
value should be increased. Additional information concerning this flowfield 
correction can be found in Appendix C. 

The artificial impingement surface is generated in the computer program 
after the potential flow calculation. Thus, the flow field is not influenced 
by the creation of the pseudo-surface. Also, this surface is discarded after 
the collection efficiency has been calculated, and is not used when a new 
ice surface is formed. 

LEQM 

This parameter is the flag to specify the initial velocity of the particle. 
In cases where an ice accretion is to be formed, the particles should be 
released in equilibrium with the flow, i.e., LEQM = 1. The option exists 
to specify the x- and y- components of the initial particle velocity when 
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LEQM = 0. This option could be used along with the options to specify 
the particle release position to simulate a droplet being ejected from a spray 
nozzle. 

LSYM 

If a body and flow field are symmetrical, the local collection efficiency 
distribution will also be symmetrical. Therefore, particle trajectory and 
impingement locations need to be calculated only for either the upper or 
lower surface. The local collection efficiency distribution is then assumed 
to be identical for the opposite surface. When LSYM = 1, the droplet im- 
pingement characteristics will be calculated only on the upper surface, and 
the local collection efficiency distribution is specified to be identical on the 
lower surface. When LSYM = 0, the droplet impingement characteristics 
are calculated for both the upper and lower surfaces. 

While some clean geometries are exactly symmetrical, the ice shapes 
rarely have exactly symmetrical surface coordinates. Forcing symmetrical 
collection efficiency distributions onto these surfaces has caused inaccurate 
ice shapes to form. Therefore, it is suggested that, unless three or fewer 
time-steps are to be performed, LSYM be set equal to 0 even for symmet- 
rical bodies at zero angle of attack. Of course, the computational time will 
be longer, but, in general, fewer problems will be encountered. 

LYOR, LXOR 

These are the x- and y-coordinate particle release flags used to indicate 
whether the particle release position will be specified by the user or de- 
termined in the computer program using the criteria discussed in Section 
4.2.6. When using the code to predict ice accretions, it is better to let the 
code determine the particle release positions. The positions that might be 
specified by the user for the clean geometry may be unsatisfactory as an 
ice accretion grows. These options have been included so that the code can 
also be used to calculate individual trajectories of particles released from a 
specific person. 
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NPL 

This parameter is used to specify the number of particle trajectories 
calculated in subroutine COLLEC to define the y 0 vs. s curve. The curve 
will be better defined when more trajectories and impingement locations 
are calculated. This, of course, is done at the expense of computational 
time. While a maximum of 50 trajectories can be calculated, NPL = 15 is 
normally specified. 

If NPL is set equal to 1, the impingement limits will not be calculated 
and the program will calculate the trajectory of only one particle. The 
particle will be released from a position specified by XORC, YORC. 

NSEAR 

The criteria to identify an impingement limit are specified by the pa- 
rameter YOLIM. If this specified parameter is specified too small, an exces- 
sive number of trajectories could be calculated. This parameter limits the 
number of trajectories that can be calculated while searching for the im- 
pingement limits; a value of NSEAR = 50 is normally input. Calculations 
of excessive numbers of trajectories can also be caused by erroneous input 
values and coordinates defining the body geometry. 

TIMSTP 

As mentioned in the description of the GEPS parameter, the integra- 
tion step size is determined in the program so that the single step error 
estimate is less than GEPS. Since the program automatically determines 
the step size during the integration, the value of the initial step size is not 
critical. A value of TIMSTP = 0.0005 has been used consistently for all 
calculations. This value is then decreased or increased as required in the 
computer program. 

NSI 

This parameter specifies the number of droplet size increments used to 
characterize the cloud droplet size distribution. If a monodispersed cloud 
is to be input, NSI = 1. The effect of using a multidispersed droplet 
distribution as opposed to a monodispersed distribution to determine an 
ice accretion shape is discussed in Section 7.4. 
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5. 2. 2. 2 Trajectory Input II (NAMELIST TRAJ2) 


CHORD 

The chord of the airfoil (or diameter for a cylindrical body) is used 
as the reference length for the coordinate release inputs discussed in this 
section. The chord is input in meters. 

G 

The acceleration of gravity is input in m/s*. If it is input as zero, 
the effect of gravity on the particle trajectories is neglected. The effect 
of gravity on the trajectories of droplets less than 50.0 microns is usually 
negligible, and is therefore omitted in most icing studies 19 . The effect of 
gravity was omitted for all of the sample cases in this report except for the 
ones using a multidispersed droplet distribution containing droplets larger 
than 50.0 microns. 

PIT, PRATK 

PIT is the initial particle pitch angle, which is defined as the angle be- 
tween the axis oriented parallel to the airfoil x- axis and the flight reference 
line (See Figure 4.3 and Section 4.2.2). For spherical, non-rotating particles 
(such as water droplets) PIT = 0.0. PRATK is the initial particle angle of 
attack and should also be set equal to 0.0 for spherical particles. 

When LEWICE is used to calculate ice accretion shapes, PIT and 
PRATK should both equal 0.0. The option to calculate the trajectories 
for non-spherical, rotating particles has been included so that the particle 
trajectory calculation may be useful for alternate applications, however, 
appropriate equations for the lift, drag, and moment coefficients must be 
supplied by the user and placed in subroutine COEFF. The subroutine was 
developed so that the coefficients are functions of the Reynolds number and 
the particle angle of attack. 

XORC, YORC 

These parameters are used to specify the particle release location when 
LXOR = 0 and LYOR = 0. XORC and YORC should be input normalized 
with respect to the airfoil chord length. 

XSTOP 


51 



This parameter is the maximum downstream value of x/chord for which 
particle trajectories are calculated. If a particle reaches a location where 
x > XSTOP, it is considered to have missed the body and moved out of 
range. This rear boundary of the computational box should extend at least 
past the location of maximum thickness, and often further, depending on 
the geometry and angle of attack. If the value of XSTOP is greater than the 
maximum x-coordinate defining the body (XREAR), XSTOP is set equal 
to XREAR. 

YOMAX, YOMIN 

These are the initial guesses for the y-coordinate of the upper and lower 
surface tangent trajectory release points, normalized with respect to the 

chord. 

YOLIM 

YOLIM is the accuracy criteria to be used in determining when an im- 
pingement limit has been reached. As discussed in Section 4.2.6, when the 
release points of two trajectories (one that hit the body and one that missed 
the body) are within YOLIM, the trajectory that hit the body is identified 
as either the upper or lower surface tangent trajectory. The smaller the 
value of YOLIM, the greater the number of trajectories that will have to be 
calculated to identify the tangent trajectory. A practical value for YOLIM 
is 0.0001; this value is used for all sample calculations. 

VXPIN, VYPIN 

These values are the x- and y-components of the initial particle velocity 
in m/s. These values need to be input only when LEQM = 0. 

5. 2.2.3 Droplet Distribution Input (NAMELIST DIST) 

DPD 

This array is used to specify the droplet sizes characterizing the distri- 
bution. The droplet sizes should be input in microns. 
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FLWC 


FLWC is the fraction of the total liquid water content contained in each 
droplet size increment. 

CFP 

This parameter is the Cunningham correction factor (Reference 15). 
Small particles, i.e., those less than 10.0 microns, have a drag slightly less 
than that given by the equation for spheres in cross-flow. The Cunning- 
ham correction factor is used to correct the drag coefficient given by the 
equations in subroutine COEFF. The correction factors are found in Table 
4.2. 


5.2.3 Ice Accretion Input (NAMELIST ICE) 

VINF, TAMB, PAMB, LWC, DPMM, RH 

These variables are used to specify the icing condition. The pressure and 
temperature inputs are static conditions. The variable DPMM is the mass 
median droplet diameter of the specified droplet distribution in microns. 
If a monodispersed cloud is specified (NSI = l), DPMM is equal to the 
droplet size specified as DPD. When a multidispersed cloud is specified 
(NSI = l), DPMM should still be input because it is used as a label on the 
parameter plots. 

XKINIT 

This parameter is the initial value of the equivalent sand-grain rough- 
ness of the icing surface. The integral boundary layer method applied to 
calculate the convective heat transfer coefficient uses this variable to ac- 
count for the effect of surface roughness. The value of XKINIT is input 
by the user and obtained using a relationship that expresses XKINIT as a 
function of static temperature, velocity, and LWC. The value of XKINIT 
is constant throughout the icing encounter. A discussion of the procedure 
used to determine this relationship is given in Appendix F. 

SEGTOL 

To maintain adequate definition of the body geometry, it is necessary 
to divide segments as they grow. This variable is the maximum fractional 
amount by which a segment may increase in length before being divided in 
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two. A value of SEGTOL < 2.0 will result in progressively shorter segments. 
Values greater than 2.0 will allow the segments to grow progressively longer. 
A value of SEGTOL = 1.5 has been found to work well in most icing 
predictions. 

QCOND 

This variable is used to input heat flow either to or from the body 
surface. The values of QCOND are input in W/m 2 as a function of surface 
distance. It can be used to simulate thermal anti-icing systems, but is not 
applicable to de-icing systems since the thermodynamics of the ice-to-liquid 
phase change at the surface of the body are not modeled. The results can 
only be assumed to be correct for the first timestep because, the effect 
of conduction through an ice layer formed during the first timestep is not 
modeled. 

5.3 Interactive Input 

The timestepping feature of LEWICE can, at times, cause conditions to 
exist that require the user’s interaction. This section describes the pri- 
mary interactive prompts and responses required in the computer code. 
Several interactive prompts are also made to indicate an unusual condition 
requiring analysis by the user. These will be discussed as they occur in the 
examples in Section 7.0. 

After the primary input file (unit 35) has been read and the potential 
flow and particle trajectory files have been set up, the program will prompt 
the user to enter the desired icing time (in seconds). The icing time is 
considered to be the total length of time that the accretion is to be grown. 
The prompt will read as follows: 

TIME = 0 

ENTER DESIRED ICING TIME (SEC), (F10.0) 

Upon entering the icing time, the program will then prompt the user to 
enter the desired time increment for the first timestep with the following 
statement: 

ENTER DESIRED TIME INCREMENT (SEC), (F10.0) 



A general guideline for selecting an icing time increment for a given icing 
condition is discussed in Section 5.3.1. 

The user will then be asked to select the desired plot options. The 
prompt will be as follows: 

AVAILABLE PLOT OPTIONS 

0 - NO PLOTS 

1 - PARAMETER PLOTS ONLY 

2 - TRAJECTORY PLOTS ONLY 

3 - PARAMETER AND TRAJECTORY PLOTS 
ENTER PLOT OPTION ( ) 

If plot option 2 or 3 is chosen, the trajectories will be plotted immediately 
after each is calculated. The trajectory points are not saved and, therefore, 
are lost when subsequent trajectories are calculated. Plot options 1 and 3 
will send the program into a plotting routine after the completion of the 
timestep so that the significant parameters can be plotted. The plotting 
routine uses a menu from which the following plots can be selected: 

01 - Iced airfoil 

02 - Particle release point (yo) vs. Surface impact distance(s) 

03 - Local collection efficieny (/?) vs. Surface distance(s) 

04 - Edge velocity (V e ) vs. Surface distance(s) 

05 - Edge temperature (T e ) vs. Surface distance(s) 

06 - Edge pressure (P e ) vs. Surface distance (s) 

07 - Surface temperature (T, ttr ) vs. Surface distance(s) 

08 - Convective heat transfer coefficient (h e ) vs. Surface distance(s) 

09 - Equivalent sand-grain roughness height (k,) vs. Surface distance(s) 

10 - Ice density (pi) vs. surface distance(s) 

11 - Freezing fraction (f) vs. surface distance(s) 

Examples of these plots are given in Section 6.0, LEWICE Output. 

5.3.1 Selection of an Icing Time Increment 

The timestepping procedure makes LEWICE unique among ice accretion 
prediction methods because it allows the physics of the ice accretion to 
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be more accurately modeled. Unfortunately, the procedure also adds com- 
plexity to the model because a proper time increment must be selected by 
the user. Also, the timestep will influence the ice accretion shape for a 
specific icing condition; the shape is not uniquely determined by the icing 
conditions. For example, Figure 5.3 shows three glaze ice accretion shapes 
calculated for the same icing condition, but with different timesteps. In 
Figure 5.3a, the accretion was formed in a single timestep. This shape 
has the basic shape of the experimental ice accretion but lacks much of 
the detail. Figures 5.3b and c show the same accretion formed at shorter 
timesteps. Note that, as the timestep is decreased, the predicted accretion 
takes on more of the characteristics of the experimental ice shape. Similar 
results for a rime ice accretion are shown in Figure 5.4. These results indi- 
cate that there is some maximum amount of ice that should be deposited 
during a single timestep. On the other hand, when a short timestep is se- 
lected, more steps are required to form the final ice shape, which increases 
the computational time required for each icing condition. 

Many ice accretion shapes have been calculated during the development 
of LEWICE, and a criterion has been developed to help the user select an 
appropriate timestep. Various authors have developed a term known as the 
accumulation parameter, which is given by the following equation: 


A c = 


fi\maz Voo LWC At 


C Pi 


(25) 


where c is the chord length of the body geometry (Some authors omit the 
factor P\maz in the definition of A™.) The accumulation parameter is rep- 
resentative of the non-dimensional maximum thickness of the ice accreted 
during time At. Therefore, a limiting value of A c could be used to de- 
termine the length of a timestep. However, it was found that when ice 
shapes were formed on an airfoil at an angle of attack, better results were 
obtained using timesteps shorter than those used when the airfoil was at 
0.0 . Equation (25) was modified as follows: 


A e = 


P\ma* V'oo (LWC) A t (1 + *s) 


C pi 


(26) 


where a is in degrees. When comparisons with experimental ice accretions, 
the modified accumulation parameter was generally less than .1. Therefore, 
the following equation can be used to calculate the time step: 
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(.1) (9.17 x 10* A) ‘ 
- fiU*. V„ LWC (1 + ft) 


( 27 ) 


Suppose that an ice accretion is to be formed on an airfoil with a chord 
of 0.3 m at the following icing condition: 


Velocity = 80.0 m/s 

LWC = 1.2 g/m 8 

Angle of attack =4.0 
Icing time =5.0 min 


If the value of /3\ maz , the maximum value of the local collection efficiency, 
is known, it should be used to evaluate Equation (27). If not, a value of 
0.80 is a reasonable upper limit and can be substituted into Equation (27). 
The equation is evaluated as follows: 


At < (.1) (9.17 x 10 Vm 8 ) (.3m) 

_ 0.80 (80m/s) (1.2jr/m s ) (1 + 53) 

< 29.85 seconds ~ 30.0 seconds 

This indicates that the initial timesteps should be no greater than 30.0 
seconds. Since Equation (27) was developed only to provide guidance in 
selecting a timestep, it is always appropriate to round the calculated time 
to a whole number. 

Depending on the icing condition and size of the geometry, a time in- 
crement may be calculated that is greater than the total icing time. In this 
case, it is recommended that the total icing time be divided into at least 
two, and perhaps three, timesteps even though Equation (27) indicated the 
ice could be accreted in one step. If a single step were used, none of the 
time dependent features of the accretion would be calculated by LEWICE. 
It was found that more realistic ice accretions are predicted if at least two 
timesteps are used. 

There is one additional note concerning the selection of a timestep. 
Depending on the icing condition, the predicted shape can become very 
convoluted and irregular after several timesteps and computational inac- 
curacies, such as multiple calculated stagnation points, may occur. If the 
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multiple effects of these approximations cause the abnormal termination 
of the code or the formation of a non-physical ice accretion, increase the 
timestep and re-run the condition. By doing so, the user can get an idea 
of the general size and shape of the accretion. By comparing this result 
to the accretion obtained using the shorter timesteps, the accuracy of the 
prediction can be evaluated. 
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Figure 5.3: Effect of icing time increment on a calcualted glaze ice accretion. 
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Chapter 6 
LEWICE Output 

The standard output from LEWICE consists of both printed and plot- 
ted output. The printed output is made available primarily for diagnostic 
purposes and examination of the calculated values. Most comparisons of 
calculated variables and ice shapes are made with the graphics routines 
provided in the code. The output discussed in this section is the output for 
the sample case presented in Section 7.1. 

6.1 Printed Output 

Printed output is produced by the potential flow, particle trajectory, and 
ice accretion portions of the computer code. All output is printed on unit 
56. The following sections describe the printed output from each of these 

portions. 

6.1.1 Printed Output from the Potential Flow Code (S24Y) 

All write statements contained in the original version of S24Y are included 
in LEWICE, however, many of these write statments have been commented 
out to reduce the amount of printed output. 

The initial output from the potential flow code is a listing of the input 
geometry coordinates, as shown in Figure 6.1a. Following the coordinate 
listing are the calculated non-dimensional surface velocities and the body 
pressure coefficients, as shown in Figure 6.1b. A description of each of the 
parameters printed on this output page is as follows: 
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Angle of attack (degress) specified by the user 
Calculated angle of zero lift (degrees) 
Calculated lift coefficient. The reference length 
is that specified by CHORD 
Chord length (meters) specified by the user to 
use as the reference length in the calculation 
of CL 

NO. OF ELEMENTS, Number of segments specified to define 
TOTAL ELEMENTS the body geometry 
I, J Segment number 

X, Y x,y coordinates of the midpoint of each 

segment 

S Surface distance to the midpoint of each 

segment 

s = 0.0 corresponds to the trailing edge of the 
body (point number 1) 

VT Non-dimensional surface velocity given by the 

following equation: 


ALPHA 
ALPHA O 
CL 

CHORD 


CP 



Surface pressure coefficient, c p calculated 
from VT using the following equation: 


CP = 1.0 - VT 2 

The force coefficients (lift, moment, normal, axial, and pressure drag) 
are calculated from the pressure coefficient data and printed after the data 
described above. These variables are calculated using two methods from 
the same set of pressure coefficient data. The first method integrates the 
pressure coefficient curve using the trapezoidal rule while the second uses 
Simpson’s rule. 
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6.1.2 Printed Output from the Particle Trajectory Code 

The first output from the particle trajectory code will be a statement iden- 
tifying the x location from which the particles were released, XO. This 
output, shown in Figure 6.2, is a result of the calculations performed in 
subroutine RELEAS. 

The next output consists of the identification of some of the geometry 
characteristics. All of the values are dimensional with standard units in the 
MKS system. Recall that the geometry coordinates have been adjusted to 
avoid the discretization errors in the flow field calculation (Appendix C). 
The values of the leading edge, trailing edge, and thickness are determined 
from these modified values and not from the original input coordinates. 
This computational correction is removed before the collection efficiency 
calculation and the calculation of the new ice surface. 

The particle trajectory data begin after the geometry characteristics are 
printed. The first output is a statement identifying how the particles were 
released. If the particles are released in equilibrium with the air (LEQM = 
1), the following message will be printed: 

The particles are released in equilibrium with the air. 

If LEQM = 0, no message is printed, and the particle release velocities 
are specified by the user. In this case, the initial particle velocities are 
printed in the line of data shown in Figure 6.2. 

Note that when LEQM = 1, the initial particle velocities are shown to 
be 0.0 m/s. This indicates that the user did not specify the input velocity. 
The initial velocity of each particle is determined by the program and will 
depend on the particle release location. 

The next statement will indicate the percent mass of each particle cor- 
responding to the droplet diameter specified. Following this statement, the 
results of the integration of the particle trajectory equations are printed. 
This output indicates whether a particle released from a point XO, Y0 im- 
pinges on the body. The column headings used in this section of the output 
are defined as follows: 



XO = x-location of particle release 

YO = y-location of particle release 

XP = x-location where particle either impinged 
upon the body or moved out of range 

YP = y-location where the particle either impinged 
upon the body or moved out of range 

S = surface distance from the stagnation point 

to the particle impingement point (lower surface 
is negative) 

DT = size of the integration timestep when the particle 
either hit the body or moved out of range 

NSTP = number of integration timesteps required for the 

particle to either impinge upon the body or move out 
of range 


As discussed in Section 4.2.6, particle trajectories are first calculated 
to determine the impingement limits in subroutine IMPLIM. The release 
points for the upper and lower surface tangent trajectories are defined as 
the particle release positions. A particle released between these points 
will strike the body, and one released outside of these points will miss the 
body. When the impingement limits are found, they will be printed as 
shown in Figure 6.2. The remaining particle trajectory calculations, for 
particles released between the upper and lower surface tangent trajectory 
release points, are made in COLLEC. All of these particles should strike 
the body, and are used to determine the local collection efficiency. When 
a droplet distribution has been specified, the output shown in Figure 6.2 
will be repeated for each droplet size before continuing with the collection 
efficiency calculation. 

After completing the calculation of the impingement limits for each 
droplet size, the program enters subroutine EFFICY, where the local col- 
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lection efficiency is calculated. The output consists of the calculated y a vs. 
s points that form the curve that is differentiated to determine the local 
collection efficiency, as shown in Figure 6.3a. This is followed by the surface 
distance and calculated local collection efficiency for each body segment. 
An example of this output is shown in Figure 6.4b. If a droplet distribution 
has been specified, the output shown in Figure 6.3a and b is printed for 
each droplet size. 

When a droplet distribution has been specified, the local collection effi- 
ciencies for each droplet size in the distribution must be combined to form 
a cumulative local collection efficiency distribution. The local collection 
efficiency for each body segment corresponding to a specific droplet size 
is weighted using the fraction of the total mass specified on input. This 
weighting procedure is described in Section 4.2.5.!. These cumulative val- 
ues are used in the thermodynamic and ice accretion portions of the code. 


6.1.3 Printed Output from the Ice Accretion Code 

The first page of output from the thermodynamic and ice accretion por- 
tions is shown in Figure 6.4a. This output contains the run identification 
specified by the user on input, and the current icing time in seconds. The 
free stream icing conditions are then printed, followed by general informa- 
tion concerning the thermodynamic and ice accretion calculations. This 
information contains the body segment numbers corresponding to the stag- 
nation point, the upper and lower surface boundary layer transition points 
(transition from laminar to turbulent flow), and the upper and lower sur- 
face icing limits. The total number of points used in the calculation of 
this timestep is given, followed by the number of segments added to the 
geometry generated by the previous timestep. 

Following this page are three pages containing the detailed output of the 
most significant aerodynamic, thermodynamic, and ice accretion parame- 
ters. All of these parameters are functions of the surface distance, s, and 
are shown in Figure 6.4b. The column headings in the computer output 
are defined as follows: 

Page 1 

I - Body segment number 

X - x-coordinate of the iced surface corresponding 


66 



Y 


to segment I (m) 

- y-coordinate of the iced surface corresponding 
to segment I (m) 

S - Surface distance, s, to the midpoint of segment I (m) 

VE - Velocity at the outer edge of the boundary 

layer, V;, (m/s) 

TE - Static temperature at the outer edge of the 

boundary layer, T„ (K) 

PE - Static pressure at the outer edge of the 

boundary layer, P e , (Pa) 

RA - Density of the air at the outer edge of the 

boundary layer, p a , (kg/m 3 ) 

SEGLENGTH - Length of the body segment, s, (m) 


Page 2 


HTC 

XK 

BETA 

FFRAC 

RI 

TSURF 

DICE 


- Convective heat transfer coefficient, h c , 
(W/m 3 /K) 

- Equivalent sand-grain roughness height, 

M 

- Local collection efficiency, 0 

- Freezing fraction, f 

- Density of the ice, pi, (kg /m 3 ) 

- Equilibrium surface temperature, T tur , (K) 

- Thickness of the ice accreted in the current 

timestep, d,, (m) 


Page 3 


QCOND 

MDOTC 

MDOTRI 

MDOTE 


- Conductive heat flux from the body surface, 
q e ,(W/m*) 

- Mass flux of liquid impinging in segment I, 

m e (kg/s) 

- Mass flux of liquid water running along the 
surface into the control volume, 

of segment I 

- Mass flux of water vapor evaporating from 
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the control volume, m e , of 
segment I (kg/s) 

MDOTTI - Total mass flux of water entering control 
volume, rhr in , of segment I (kg/s) 

MDOTT - Total mass flux of water in to control volume 

ririT, of segment I (kg/s) 


This concludes the description of the printed output from LEWICE. All 
of the output described above is repeated for each time step. 


6.2 Graphical Output 

While the printed output is useful to verify whether the code is performing 
as expected, much of the output from LEWICE is displayed graphically to 
aid in the evaluation of the results. 

The plotting commands used in LEWICE are unique to NASA Lewis 
Research Center and, therefore, are not likely to be directly applicable to 
another facility. GRAPH2D/GRAPH3D commands are used. 

All plots discussed in this section have the same border , which contains 
title and icing condition information. 

6.2.1 Droplet Trajectory Plots 

As discussed in Section 5.3.1, when plot option 2 or 3 is selected, the particle 
trajectories will be plotted. An example of five such plots are shown in 
Figure 6.5. The axes limits (in meters) are determined by the program 
so that the entire body geometry is plotted. When a particle impinges 
upon the body, the y-coordinate of the release point, yo, and the surface 
impingement distance, s, are displayed on the plot. If the particle misses 
the body, this information is omitted. 

The plotting commands are located in subroutine PLTRAJ, which is 
called from subroutine INTIG after the calculation of the trajectory is com- 
pleted. Each trajectory is plotted immediately after it is calculated and, 
once the calculation of a new trajectory is begun, the previous trajectory 
coordinates are erased. 
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6.2.2 Ice Accretion Data 


Much of the ice accretion data (Section 6.1.3) can also be plotted to as- 
sist the user in interpreting the results. If plot option 1 or 2 has been 
selected, the program will enter the plotting routine (subroutine PLOTD) 
after completing the ice accretion calculations. The following plot menu 
will be displayed upon entering the plot routine: 

AVAILABLE PLOT OPTIONS 

0 - NO PLOTS 

1 - ICED GEOMETRY 

2 - Z VS S 

3 - BETA VS S 

4 - VE VS S 

5 - TE VS S 

6 - PE VS S 

7 - TSURF VS S 

8 - HTCVSS 

9 - XK VS S 

10 - ICE DENSITY VS S 

11 - FFRAC VS S 
ENTER OPTION NUMBER (12) 


The user then inputs the two-digit identifier for the desired plot. 

When plot 01 is selected the iced geometry, with all preceding timesteps, 
will be plotted. Before plotting the geometry, the user will be asked to 
specify the percent of the geometry to be plotted with the following prompt: 

ENTER PERCENT OF GEOMETRY TO BE PLOTTED (F10.0) 

Since the size of the plot is fixed, the larger the portion of the airfoil to 
be plotted, the smaller the geometry will appear on the plot. For example. 
Figures 6.6a, b, and c show geometries plotted with the specified percent 
equal to 100., 50., and 25., respectively. After the geometry is plotted, the 
program will return to the plot menu. 

Plots 02 to 11 are all parameters that are functions of the surface dis- 
tance, s. The format for all of these plots is similar and is described below. 
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After selecting the desired plot, the minimum and maximum values on 
the ordinate and abcissa will be displayed. After each maximum and min- 
imum is displayed, the user will be asked to specify the desired axis limit. 
After the desired limits are input, the program will ask if experimental data 
is to be plotted, i.e., 

ENTER 1 IF EXPERIMENTAL DATA IS TO BE PLOTTED IF NOT, 
ENTER O 

If experimental data are to be plotted, the program will first ask for the 
number of data points, and then ask the user to manually input the data. 
After the last experimental data point is input, the plot will be displayed, 
and the experimental data will be represented as circles. 

After viewing the plot, the user will be given the chance to change the 
axes limits with the following prompt: 

IF YOU WOULD LIKE A DIFFERENT SCALE, ENTER 1. IF NOT, 
ENTER O. 

If this is desired, the maximum and minimum ordinate and abcissa values 
will be displayed again, and the above procedure is repeated. If new axes 
limits are not desired, the plot menu will be returned to the screen. 

As previously stated, plots 02 to 11 all use the operating format de- 
scribed above with the occasional exception of plot 02, y 0 vs. s. When a 
droplet distribution is specified, the y 0 vs. s curve corresponding to each 
droplet size will be plotted. The correct droplet size for each plot will be 
printed in the parameter box in the upper left-hand corner of the plot. The 
user will be asked to input the maximum and minimum axes limits for the 
plot of each droplet size. 

Examples of each of the parameter plots (02-11) are shown in Figures 
6.7a-k. 
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UNTP. AMS FORMED 

COORDINATE DATA 

FOR BODY ID = 

1, NACA 

0012 : EXAMPLE 

x 

X(I) 

Y(I) 

I 

XU) 

Y(I) 

1 

1 .0000000 

0.0000000 

51 

0.0150000 

-0 . 0207970 

2 

0 . 9899999 

-0.0018740 

52 

0.0100000 

-0 .0172480 

3 

0 . 9800000 

-0.0036110 

53 

0.0075000 

-0.0150780 

4 

0.9700000 

-0.0052390 

54 

0.0050000 

-0 .0123860 

5 

0 . 9500000 

-0 . 0082510 

55 

0.0037500 

-0 .0106990 

6 

0 . 9250000 

-0.0117000 

56 

0.0025000 

-0 .0086200 

7 

0.9000000 

-0.0149080 

57 

0.0022500 

-0.0081360 

ft 

0 .8750000 

-0.0179530 

58 

0.0020000 

-0.0076230 

9 

0.8500000 

-0 . 0208880 

59 

0.0017500 

-0 .0070750 

1 0 

0 .8250000 

-0.0237390 

60 

0.0015000 

-0.0064850 

1 1 

0 .8000000 

-0 . 0265150 

61 

0.0012500 

-0.0058470 

1 7 

0 . 7750000 

-0.0292210 

62 

0.0010000 

-0.0051460 

i 3 

0 . 7500000 

-0.0318550 

63 

0.0008750 

-0.0047660 

14 

0 . 7250000 

-0.0344110 

64 

0.0007500 

-0 .0043630 

1 5 

0 . 7000000 

-0.0368870 

65 

0.0006250 

-0 .0039310 

1 s 

3 . 6 750 0 00 

-0 .0392810 

66 

0.0005000 

-0 .0034620 


3 . 6 5 0 0000 

-0 .0415850 

67 

0.0003750 

-0.0029450 

18 

0 .6 250 00 0 

-0 . 0437950 

68 

0.0002500 

-0.0023560 

1 9 

0 .600 0 000 

-0 . 0459040 

69 

0.0001250 

-0.0016320 

20 

0 . 57 50 000 

-0 . 0479060 

70 

0.0000000 

0 . 0000000 

21 

0 . 5500000 

-0 .0497930 

71 

0.0001250 

0 . 0016320 

22 

0 . 5250 0 00 

-0.0515600 

72 

0.0002500 

0 .0023560 

23 

0.5000000 

-0 .0531980 

73 

0.0003750 

0 .0029450 

24 

0.4750000 

-0 . 0546980 

74 

0.0005000 

0 .0034620 

25 

0.4500000 

-0.0560510 

75 

0.0006250 

0.0039310 

26 

0 . 4250000 

-0.0572430 

76 

0.0007500 

0.0043630 

27 

0 .4000000 

-0 . 0582620 

77 

0.0008750 

0 .0047660 

28 

0 . 37 50000 

-0 .0590900 

78 

0.0010000 

0.0051460 

29 

0 . 3500000 

-0.0597070 

79 

0.0012500 

0.0058470 

30 

0 . 3250000 

-0.0600940 

80 

0.0013000 

0.0064850 

31 

0 . 3000000 

-0.0602260 

81 

0.0017500 

0.0070750 

32 

0.2750000 

-0.0600760 

82 

0.0020000 

0.0076230 

J3 

0.2500000 

-0.0596120 

83 

0.0022500 

0.0081360 

34 

0.2250000 

-0.0587940 

84 

0.0025000 

0.0086200 

35 

0 .2000000 

-0.0575700 

85 

0.0037500 

0 .0106990 

36 

0 .1 7 50000 

-0.0558790 

86 

0.0050000 

0.0123860 

37 

0 .1 500000 

-0.0536360 

87 

0.0075000 

0.0150780 

38 

0 .1 250000 

-0.0507320 

88 

0.0100000 

0.0172480 

39 

0 .1 00 0000 

-0.0470040 

89 

0.0150000 

0.0207970 

40 

0 . 090000 0 

-0.0452280 

90 

0.0200000 

0 .0237260 

41 

0 . 0800 000 

-0.0432520 

91 

0.0250000 

0.0262300 

42 

0 .0700000 

-0.0410380 

92 

0.0300000 

0.0284620 

4 3 

0 . 06 00000 

-0.0385350 

93 

0.0350000 

0.0304960 

44 

0.0500000 

-0.0356740 

94 

0.0400000 

0.0323620 

45 

0 . 0450000 

-0.0340520 

93 

0.0450000 

0.0340820 

46 

0.0400000 

-0 .0323620 

96 

0.0500000 

0 . 0356740 

4 7 

0.0350000 

-0 .0304960 

97 

0.0600000 

0.0385350 

48 

0 .0 30000 0 

-0 . 0284620 

98 

0.0700000 

0.0410380 

49 

0 . 025000 0 

-0.0262300 

99 

0.0800000 

0.0432520 

50 

0 . 020 000 0 

-0 . 0237260 

100 

0.0900000 

0.0452280 


I 


XCI) 

101 

0 

.1000000 

102 

0 

. 1 250 00 0 

103 

0 

. 1 500 0 0 0 

104 

0 

.1750000 

105 

0 

. 200000 0 

106 

0 

.2250000 

107 

0 

.2500000 

103 

0 

.2750000 

109 

0 

.3000000 

no 

0 

. 3250 0 0 0 

in 

0 

. 350 00 0 0 

112 

0 

. 3750 0 0 0 

113 

0 

.4000000 

114 

0 

.4250000 

115 

0 

.4500000 

116 

0 

.4750000 

117 

0 

.5000000 

118 

0 

. 52500 0 0 

119 

0 

. 55000 00 

120 

0 

.5750000 

121 

0 

.6000000 

122 

0 

.6250000 

123 

0 

.6500000 

124 

0 , 

.6750000 

125 

0 , 

.7000000 

126 

0 , 

. 7250000 

127 

0 , 

. 750 000 0 

128 

0 , 

, 7 7 50000 

129 

0 . 

.8000000 

130 

0 . 

,8250000 

131 

0 . 

,8500000 

1 32 

0 . 

,8750000 

133 

0 . 

, 9000000 

134 

0 . 

9250000 

135 

0 . 

9500000 

136 

0 . 

9700000 

137 

0 . 

9500000 

138 

0 . 

9899999 

139 

1 . 

0000000 


a. Geomentry coordinates input. 


ORIGINAL 
OF POOR 


DOUGLAS AIRCRAFT COMPANY TWO-DIMENSIONAL POTENTIAL FLOW PROGRAM 


COMBINED 

FLC14 

NACA 0012 : EXAMPLE 1 

ALPHA = 

0.000000 

ALPHA 0 = 

-0.000010 

CL = 

-0.000001 

CHORD 

1.000000 

BODY ID ■ 

= 1 

NACA 0012 : EXAMPLE 1 

I 

X 

Y 

s 

1 

0 . 9950 0 30 

-0 .000 9536 

0 . 0024939 

2 

0 . 9850 022 

-0.0027573 

0.0074755 

3 

0 .97 500 18 

-0 .0 044 365 

0 .0124468 

4 

0.9600045 

-0 .0067 780 

0 .0 198878 

5 

0 . 9 37 5044 

-0 . 0100112 

0.0310313 

6 

0 .9 1 250 30 

-0.0133284 

0 . 0433954 

7 

0 .8875018 

-0 .0164470 

0 . 0557466 

5 

0.8625011 

-0 ,0 1 94323 

0 . 0680897 

9 

0 .8 3 7 500 7 

-0 .0223233 

0 . 0804271 

10 

0 .8 1 2500 5 

-0 .0251 359 

0.0927601 

1 1 

0 . 78 7 5008 

-0.0278768 

0 . 1 0508 92 

12 

0 . 7 6 25008 

-0 .0305473 

0.1174146 

13 

0.7375007 

-0 .0331427 

0 . 129736 1 

14 

0 . 7 1^5006 

-0 .0356590 

0 . 1420537 

15 

0 . 68 750 0 5 

-0 . 0 380 946 

0 .1 543674 

16 

0 .66250 0 9 

-0 . 04 Q4444 

0 . 1666770 

17 

3 . 6 3 7 5 0 08 

-0 . 042 7 021 

0 .1789825 

18 

0.6125008 

-0 . 0448624 

0.1912838 

19 

0 . 587 500 7 

-0.0469188 

0.2035809 

20 

0 . 56 25006 

-0 . 048864 1 

0.2158735 

21 

0 . 53 7 50 1 0 

-0 . 0506 920 

0.2281619 

^ -9 

(p c 

0 .5125009 

-0 . 0523956 

0.2404459 

23 

0 . 487 50 1 0 

-0.0539657 

0.2527257 

24 

0 . 46250 1 0 

-0 . 0553936 

0.2650014 

25 

0 .437 500 9 

-0 .0566678 

0.2772729 

26 

0 .4125009 

-0.0577752 

0.2895406 

27 

0 . 3875008 

-0 . 0587 0 1 0 

0 . 3018047 

28 

0 . 36 250 06 

-0 . 0594260 

0.3140656 


NO. OF BODIES 1 
TOTAL ELEMENTS 138 


NO. OF ELEMENTS 138 


VT 

CP 

J 

SIGH A 

-0.8268576 

0.3163066 

1 

-0.013142 

-0.8886778 

0.2102517 

2 

-0.013130 

-0.9189486 

0.1555335 

3 

-0.012671 

-0 . 9451062 

0 . 1 06774 3 

4 

-0 .0 1200 1 

-0.9680851 

0.0628113 

5 

-0.011208 

-0.9839670 

0.0318090 

6 

-0 . 0 1 056 0 

-0.9951780 

0 . 00 96207 

7 

-0.010141 

-1.0043350 

-0.0087891 

8 

-0 .00 9908 

-1.0128355 

-0.0258350 

9 

- 0 .0 0 97 6 1 

-1.0208483 

-0 .042 1 305 

10 

-0.009632 

-1.0285645 

-0.0579443 

1 1 

-0.009518 

-1.0361195 

-0 .0735426 

12 

- 0 . 0 0 9385 

-1 . 0434093 

-0.0887022 

1 3 

-0.009215 

-1.0504131 

-0 . 1033669 

14 

-0.009038 

-1.0573263 

-0.1179380 

15 

-0.008855 

-1.0641356 

-0.1323843 

16 

-0.008629 

-1 .0707989 

-0 .14660 93 

17 

-0 . 0 08389 

-1.0773840 

-0 .160756 1 

18 

- 0 . 0 08 1 16 

-1.0838785 

-0 . 1747923 

19 

-0 .0 078 18 

-1.0902719 

-0 . 1886921 

20 

-0.007483 

-1.0966110 

-0.2025557 

21 

-0.007134 

-1.1029654 

-0.2165318 

22 

-0.006745 

-1.1092901 

-0.2305241 

23 

-0.006321 

-1.1156693 

-0.2447176 

24 

-0.005863 

-1 . 1220798 

-0.2590628 

25 

-0.005340 

-1.1285248 

-0.2735682 

26 

-0 . 004 7 72 

-1.1350374 

-0.2883091 

27 

-0.004120 

-1.1414566 

-0.3029222 

28 
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b. Surface flow characteristics 


Figure 6.1: Printed output for the potential flow calculations in Example 
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Surface flow characteristics (continued) 


Figure 6.1: continued 
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Surface flow characteristics (continued) 


Figure 6.1: Concluded 
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Figure 6.2: Printed output for the individual particle trajectory calculations 
in Example 1. 
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a. Calculate y 0 vs s points 
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b. Local collection efficiency for each body segment. 


Figure 6.3: Summary of the particle trajectory and local collection efficiency 
calculations in Example 1. 
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b. Local collection efficiency for each body segment (continued). 


Figure 6.3: Concluded 
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a. Icing condition and summary of ice accretion data 


Figure 6.4: Printed output from the ice accretion calculations in Example 
1. 75 
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b. Calculated surface characteristics (continued) 
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NACA 0012 : EXAMPLE 1 

ICING CONDITION: 


: TIME = 60.000 SEC 


STATIC TEMPERATURE <C) 260.55 
STATIC PRESSURE (PA) 90748.00 
VELOCITY (M/S) 129.00 
LUC (G/M**3) 0.50 
DROPLET DIAMETER (MICRONS) 20.00 
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b. Calculated surface characteristics (continued) 


Figure 6.4: continued 
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b. Calculated surface characteristics (continued) 


Figure 6.4: continued 
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b. Calculated surface characteristics (continued) 


Figure 6.4: Concluded 
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b. Trajectory of particle released from YOMIN 
Figure 6.5: Examples of droplet trajectory plots from Example 1. 
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e. Trajectory calculated while searching for the upper impingement limit. 

Figure 6.5: Concluded 
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c. Percent plotted = 25% 
Figure 6.6: Concluded 


85 





SCIA 


ORIGINAL PAGE IS 
OF POOR QUALITY 


VELOCITY (M/S) 129.96 

TEMPERATURE <C> -12.60 

PRESSURE (KPA) 90.75 

HUMIDITY (*) 100.00 

LHC (G/M » « 3 ) 0.50 

DROP 0 I A M (MICRON) 20.00 
TIME (SEC) 60.00 




a. Iced airfoil geometry 



S <H) 


b. Particle release position vs. surface 
impact distance 



c. Local collection efficiency vs. surface 
location 


s cm: 

d. Edge velocity vs. surface location 


c r 



- ! 6 *- 

* 2 3/35 - .'3 5 O: -.V-- D. : Ji 0/59 3. J 35 

s (m: 

e. Edge temperature vs. surface location 


:io - 



9 ? .' s 5 " - .1 S '" ~ T a i ~ - ; ! i t - :z T~r r j i r~ a rr r-g ' . : i + : 3 5 

S CM) 

f. Edge pressure vs. surface location 


Figure 6.7 Icing parameter plots for the first timestep of Example 1. 
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Chapter 7 
EXAMPLE CASES 


In this section, the capabilities of LEWICE are illustrated through the 
presentation of five example cases. These examples include the input data 
sets and the printed and plotted output data; however, printed output is 
included for example 1 only. 

7.1 Example 1: Glaze Icing on a NACA 0012 Airfoil, a = 0.0° 

In this example, a 2-minute glaze ice accretion will be formed on a NACA 
0012 airfoil at an angle of attack of 0.0 degree. The accretion will be formed 
in two 1-min timesteps. The icing conditions are shown below: 


Velocity 

Static Temperature 
Static Pressure 
LWC 

Droplet Diameter 
Roughness Height 


= 129.46 m/s 
= 260.55 K 
= 90748.0 Pa 
= 0.50 g/m s 

= 20.0 microns (monodispersed cloud) 
= 0.00035 m 


The input file and the interactive input for the first timestep in Example 
1 are shown in Figures 7.1a and b, respectively. Plot option 3 was selected 
in the interactive input, therefore, the particle trajectories and icing pa- 
rameters were plotted. Figures 6.5a-e show the first five particle trajectory 
plots from this example case. Figures 6.5a and b show the trajectories of 
particles released from Y0MAX and Y0MIN. As discussed in Chapter 5, 
these trajectories are calculated to identify the release point of a particle 
that passes above the body and one that passes below the body. These 
calculations are performed in subroutine RANGE. 

After locating an upper and lower boundary trajectory, the search for 
the upper and lower impingement limit is begun by releasing a particle 
from a position halfway between the release points of the two trajectories 
identified in Figures 6.5a and b. When a particle impinges on the body, 
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the particle number (IW), particle release ordinate (YO), and the surface 
impact distance (S) are displayed, as shown in Figure 6.5c. Using the 
Newton iteration scheme described in Chapter 4.2.6, the search for the 
upper impingement limit is continued, as shown in Figures 6.5d and e. 
These calculations are performed in subroutine IMPLIM. When plot option 
3 is selected, all of the trajectories will be plotted. Only the first five have 
been shown in this example. 

After calculating the particle trajectories, the program will prompt the 
user to preview the yo vs. s data or to continue with the collection efficiency 
calculation with the following statement. 

ENTER 1 TO PREVIEW YO VS S DATA 

ENTER 0 TO CONTINUE WITH THE COLLECTION EFFICIENCY 
CALCULATION (II) 

If option 1 is selected and the axes limits are specified as follows: 

smin = -0.04 smax = 0.04 

y 0 min = -0.01 y 0 max = 0.01 


the y 0 vs. s points calculated in subroutine COLLEC will be plotted as 
shown in Figure 7.2. Upon producing the plot, the following statement will 
be printed on the screen: 

THE CALCULATED SURFACE IMPACT DISTANCES (S) AND 
RELEASE POINTS (Y0) ARE AS FOLLOWS 


I 

Y0 

S 

I 

Y0 

s 

1 

0.00850 

0.03378 

10 

- 0.00109 

- 0.00178 

2 

0.00765 

0.02078 

11 

- 0.00218 

- 0.00355 

3 

0.00656 

0.01489 

12 

- 0.00327 

- 0.00555 

4 

0.00546 

0.01090 

13 

- 0.00436 

- 0.00797 

5 

0.00437 

0.00801 

14 

- 0.00545 

- 0.01085 

6 

0.00328 

0.00557 

15 

- 0.00654 

- 0.01484 

7 

0.00219 

0.00357 

16 

- 0.00763 

- 0.02072 

8 

0.00110 

0.00180 

17 

- 0.00848 

- 0.03167 

9 

0.00001 

0.00001 

18 

0.00000 

0.00000 


HOW MANY TRAJECTORIES ARE TO 3E DELETED BEFORE 
THE COLLECTION EFFICIENCY CALCULATION? (12) 
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At this time, y 0 vs. s points can be removed, if necessary. Non- 
physical impact locations can occur when inaccurate flow fields are cal- 
culated around the horns of glaze ice accretions. In this example, it is not 
necessary to remove any points, so option 0 was selected. 

The title of the run, icing condition parameters, ice accretion data, and 
parameter plot menu are then printed on the screen as follows: 

NACA 0012 : EXAMPLE 1 : TIME = 60.000 SEC 


ICING CONDITION: 

STATIC TEMPERATURE (C) 260.55 

STATIC PRESSURE (PA) 90748.00 

VELOCITY (M/S) 129.00 

LWC (G/M 3 ) 0.50 

DROPLET DIAMETER (MICRONS) 20.00 

ICE ACCRETION DATA 

STAGNATION POINT 70 

TRANSITION POINTS (LOWER, UPPER) 68 71 

ICING LIMITS (LOWER, UPPER) 42 97 

NUMBER OF POINTS 139 

NUMBER OF SEGMENTS ADDED 0 


AVAILABLE PLOT OPTIONS 

0 - NO PLOTS 

1 - ICED GEOMETRY 

2 - Z VS S 

3 - BETA VS S 

4 - VE VS S 

5 - TE VS S 

6 - PE VS S 

7 - TSURF VS S 

8 - HTC VS S 

9 - XK VS S 

10 - ICE DENSITY VS S 

11 - FFRAC VS S 
ENTER OPTIONS NUMBER (12) 
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All of the parameter plots for the first time-step are shown in Figures 
6.7a-k. When the option for no plots is selected, the computer will respond 
with the following message: 

TIME STEP COMPLETE: TIME = 60.0 

ICING TIME INPUT HAS BEEN REACHED 

PROGRAMS OPTIONS 

1 - CONTINUE ICING, USE PREVIOUS FLOW FIELD 

2 - CONTINUE ICING, CALCULATE NEW FLOW FIELD 

3 - TERMINATE PROGRAM 


Since a second 1-min timestep with a new flow field calculation is desired, 
option 2 is selected. 

The second timestep proceeds in a manner similar to the first. The 
plots, computer prompts, and responses are shown in the order of their 
occurrence in Figure 7.3. The printed output from the first timestep in 
this example is shown in Figures 6.1-6.4. The icing parameter plots for the 
second timestep are shown in Figure 7.4a-k. The printed output for the 
second timestep is shown in Figure 7.5. 

The overall accuracy of LEWICE is checked by comparing calculated 
and experimental ice accretion shapes. The experimental ice accretion for 
the condition was obtained by Gent (Reference 20), and is compared to the 
calculated shape in Figure 7.6. 

7.2 Example 2: Glaze Icing on a NACA 0012 Airfoil, a = 4.0° 

This example is identical to Example 1 except that the airfoil is now at a 
4.0 angle of attack. The icing time for this condition is also 2-minute and 
the accretion is formed in two 1-minute timesteps. The input file is shown 
in Figure 7.7. Note that the value of CLT in namelist S24Y is now equal 
to 4.0. 

Figures 7.8 show the particle trajectories calculated in subroutine 
RANGE, where the program searches for a trajectory that passes above 
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and below the body. More accurate selections of YOMIN and YOMAX 
would have resulted in fewer trajectories being calculated while the code 
searches for the two boundary trajectories. The first three trajectories cal- 
culated in subroutine IMPLIM are shown in Figures 7.9. These trajectories 
are calculated while the program was searching for the upper surface im- 
pingement limit. No additional trajectories will be shown for this example, 
but all trajectories will be plotted when plot option 3 is selected. The y 0 
vs. s points calculated in subroutine COLLEC are shown in Figure 7.10. 
After the plot is completed, the following prompts are printed to the screen 
to allow the user to remove any incorrect points: 


THE CALCULATED SURFACE IMPACT DISTANCES (S) AND 


RELEASE 

I 

POINTS (Y0) ARE AS FOLLOWS 
Y0 S I Y0 

S 

1 

-0.11162 

0.01610 

10 -0.12296 

-0.00836 

2 

-0.11263 

0.01042 

11 -0.12425 

-0.01156 

3 

-0.11392 

0.00727 

12 -0.12554 

-0.01528 

4 

-0.11521 

0.00473 

13 -0.12683 

-0.01979 

5 

-0.11650 

0.00276 

14 -0.12812 

-0.02534 

6 

-0.11779 

0.00081 

15 -0.12942 

-0.03270 

7 

-0.11908 

-0.00103 

16 -0.13071 

-0.04386 

8 

-0.12037 

-0.00315 

17 -0.13171 

-0.06252 

9 

-0.12167 

-0.00561 

18 0.00000 

0.00000 


HOW MANY TRAJECTORIES ARE TO BE DELETED BEFORE 
THE COLLECTION EFFICIENCY CALCULATION? (12) 


As in Example 1, option 0 was selected because it was not necessary to 
remove any points. The program then entered displays the plot menu to 
the screen, and the plots shown in Figures 7.11a-i were obtained. The 
results of the first timestep are similar to those obtained in Example 1. 

After plot option 0 was selected, the prompts for the program option 
were printed, and, as in Example 1, option 2 was selected. The plot option 
was not changed from the first timestep and the calculation of the potential 
flow field was begun. When the program entered subroutine STAG, it was 
found that two “stagnation points” had been calculated. This is indicated 
by the following message and prompt: 
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2 STAGNATION POINTS HAVE BEEN CALCULATED. ENTER 1 
TO SHOW THE LOCATIONS OF THE CALCULATED STAGNATION 
POINTS. ENTER 0 TO DISPLAY OPTIONS. 

When this situation occurs, it may be necessary to form a pseudo-surface 
over the region where the multiple stagnation points exist. The procedure 
to form the pseudo-surface is described in the following section. Appendix C 
discusses additional causes for the calculation of multiple stagnation points. 

7.2.1 Application of the Pseudo-Surface 

In normal operation, select option 1 to display the locations of the stagna- 
tion points. The locations of the stagnation points are indicated by X’s, as 
shown in Figure 7.12. After displaying the plot, the x, y, and s locations 
of the stagnation points will be displayed, and the user will be asked either 
to select one of the points for use as the stagnation point in the remaining 
calculations, or to form the pseudo-surface. 

THE FOLLOWING MULTIPLE STAGNATION POINTS HAVE 
BEEN CALCULATED 

IX Y S 

57 - 0.00089 - 0.00336 0.00123 

59 - 0.00111 -0.00309 0.00172 

AVAILABLE OPTIONS 

0 - CONTINUE CALCULATIONS WITH MANUALLY-SELECTED 

STAGNATION POINT 

1 - RECOMPUTE FLOW FIELD USING A PSEUDO-SURFACE 


In this example, option 1 was selected. 

In Figure 7.12, the axes scales were such that it was difficult to identify 
the locations of the stagnation points. Therefore, the user will be asked if 
the plotting scales should be changed to enlarge the region of interest. In 
this example, the axes limits were changed as shown below: 

x t max — Umax = 0.02 

x min = tfmin ~ 0.02 
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The locations of the stagnation points will then be re-plotted using the new 
axes limits, as shown in Figure 7.13. 

The circular arc that will be placed over the multiple stagnation points 
is defined by specifying two points on the arc and the radius. The two 
points on the arc are the locations where the pseudo-surface intersects the 
actual ice surface. The best results are obtained when the ends of the 
pseudo-surface arc intersect on a tangent to the ice surface. Figure 7.13 
shows that this may be possible near y values near -0.006 and -0.001. These 
values are input by the user as shown below: 

ENTER THE Y COORDINATE OF THE DESIRED LOWER 
LIMIT OF THE PSEUDO-SURFACE (F10.0) 

-.006 

ENTER THE Y COORDINATE OF THE DESIRED UPPER 
LIMIT OF THE PSEUDO-SURFACE (F10.0) 

-.001 

The computer will respond with the two body coordinates closest to the 
desired upper and lower limits as shown below: 

THE POINTS CLOSEST TO THE LOWER LIMIT ARE 
X= 0.0019641 Y= -0.0065025 1= 52 

X= 0.0012161 Y= -0.0055978 1= 53 

THE POINTS CLOSEST TO THE UPPER LIMIT ARE 
X= -0.0028274 Y= -0.0010464 1= 74 

X= -0.0029714 Y= -0.0008570 1= 75 


The user is then prompted to enter the segment number of the desired upper 
and lower limits. In this example, segments 52 and 75 were selected as the 
lower and upper limits of the pseudo-surface, respectively. The program 
will then request that a radius be input. In this example the radius is 0.03. 
The input of these values is shown below. 

ENTER THE I VALUE OF THE DESIRED LOWER AND UPPER 
LIMITS (213) 52 75 

ENTER THE DESIRED RADIUS (F10.0) .03 

Because the selection of an appropriate radius is a manual iteration, there 
are no general criteria for selecting a radius. If an arc with the desired 


94 



radius cannot be formed between the two points specified, the program will 
respond by asking for another radius. The iced airfoil with the pseudo- 
surface will then be plotted, as shown in Figure 7.14. After the plot is 
completed, the following menu is displayed: 

THE FOLLOWING OPTIONS ARE NOW AVAILABLE (11) 

0 - SPECIFY NEW UPPER AND LOWER LIMITS 

1 - SPECIFY A NEW RADIUS 

2 - PSEUDO-SURFACE SATISFACTORY, CONTINUE 

CALCULATIONS 

3 - PSEUDO-SURFACE UNSATISFACTORY, CONTINUE 

CALCULATIONS 

WITH THE ACTUAL ICE SURFACE 

At this time, new upper and lower limits or a new radius can be selected, 
and the above procedure will be repeated. Option 2 will cause the pro- 
gram to re-calculate the flow field using the pseudo-surface. Option 3 is 
selected when a satifactory pseudo-surface cannot be formed and any errors 
in the flow field must be accepted. In this example, the pseudo-surface was 
satisfactory, therefore, option 2 was selected. 

While it did not occur in this example, multiple stagnation points may 
still be calculated after the pseudo-surface has been formed. Unless a 
smoother pseudo-surface can be formed, it is usually best to select one 
of the points and continue with the calculations. Select the point that is 
closest to the stagnation point calculated in the previous timestep. This 
point can be identified by the value of surface distance, s. (Recall that s = 
0.0 denotes the stagnation point.) 

In this example, a single stagnation point was calculated when the 
pseudo-surface was placed on the iced surface, and the second timestep 
proceeded in a manner similar to the first. It was not necessary to remove 
any y 0 vs. s points from the plot shown in Figure 7.15 when the following 
prompt was displayed: 

HOW MANY TRAJECTORIES ARE TO BE DELETED BEFORE 

THE COLLECTION EFFICIENCY CALCULATION? (12) 

The plots obtained from the plot menu for the second timestep are 
shown in Figures 7.16ari. The calculated airfoil shape is compared to the 
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experimental ice accretion shape, obtained from Reference 20, in Figure 
7.17. 


7.3 Example 3: Rime Icing on an NACA 0012 Airfoil, a — 0.0° 


In this example, ice will be accreted on an NACA 0012 airfoil at an angle 
of attack of 0.0 for 2 min. As in Examples 1 and 2, the accretion will 
be formed in two 1-min time steps. The icing conditions, which should 
produce a rime ice accretion, are as follows: 


Velocity 

Static Temperature 
Static Pressure 
LWC 

Droplet Diameter 
Roughness Height 


= 64.73 m/s 
= 260.55 K 
= 90748.0 Pa 
= 0.50 g/m s 

= 20.0 microns (monodispersed cloud) 

= 0.00025 m 


The input file for this example is shown in Figure 7.18. 

The program executes in a manner similar to Example 1. The parameter 
plots for the first and second timesteps are shown in Figures 7.19a-i and 
7.20a-i, respectively. The predicted ice accretion shape is compared to the 
experimental shape (Reference 20) in Figure 7.21. 

As discussed in Section 4.3.1, the difference between rime and glaze ice 
is determined by the calculated value of the freezing fraction, f (FFRAC in 
LEWICE). Figure 7.19k shows that the freezing fraction was equal to 1.0 
over the entire surface for the first timestep. This would be indicative of a 
solid rime ice accretion. On the second timestep, the freezing fraction was 
slightly less than 1.0 near the stagnation point, as shown in Figure 7.20k. 
This indicates that the accretion is starting to take on some glaze ice accre- 
tion characteristics in that region, and would therefore constitute a mixed 
ice accretion. Glaze characteristics become more evident as the freezing 
fraction approaches 0.0. Experimentally, these initial glaze characteristics 
would probably not be observed until the ice accretion grew larger and the 
size of the glaze portion increased. 
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7.4 Example 4: Specification of a Droplet Size Distribution 

The atmospheric conditions in this example are identical to those in Ex- 
ample 1. However, in this example, a droplet size distribution was used to 
characterize the icing cloud instead of a single droplet size. The icing time 
for this case is 60.0 sec, and the ice will be accreted in a single timestep. 

The droplet size distribution was specified to have a normal volume 
(mass) distribution with a mean of 20 microns and a = 5, as shown in Figure 
7.22. As discussed in Section 4.2.5.1, the droplet size distribution is input 
into the code by specifying the mass fraction corresponding to a discrete 
droplet diameter. The procedure for quantifying a known distribution for 
input into LEWICE is discussed below. 

The volume fraction (fraction LWC) corresponding to each droplet size 
used to produce Figure 7.22 is shown in Table 7.1. The cumulative volume 
fraction (CVF) is calculated for each droplet size, d,, using the equation 

CVFi = CVF V . „ + (29) 

The resulting plot of CVF vs. d is shown in Figure 7.23. For input 
into LEWICE, the droplet size range must be divided into a number of 
discrete droplet size bins of a constant width. The size of the bin width is 
arbitrary, but no more than 10 bins can be input to the code. Using more 
bins will increase the accuracy of the calculations but will also increase 
computational time. Five bins with a width of 8.0 microns were used in 
this example. The droplet diameter corresponding to the middle and the 
right-hand side of each bin is shown in Table 7.2. 
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Table 7.1 Volume Fraction for a Normal Mass Distribution of 
Droplet Size with a Mean of 20mm and <7 = 5 

d (microns) Volume Fraction 


0 

0.0000 

2 

0.0002 

4 

0.0010 

6 

0.0032 

8 

0.0090 

10 

0.0216 

12 

0.0444 

14 

0.0776 

16 

0.1158 

18 

0.1474 

20 

0.1596 

22 

0.1474 

24 

0.1158 

26 

0.0776 

28 

0.0444 

30 

0.0216 

32 

0.0090 

34 

0.0032 

36 

0.0010 

38 

0.0002 

40 

0.0000 


The cumulative volume fraction at the right-hand side of each bin must 
be determined using the data in Table 7.1 and Figure 7.23 using the fol- 
lowing equation: 


V). = CVFi - CVFi _i (30) 

This volume fraction is assigned to the droplet size at the middle of the bin 
and assumes a uniform distribution of droplets inside the bin. The values 
of Vf { for each droplet size bin are shown in Table 7.2. This data is input 
into LEWICE using namelist DIST, as shown in Figure 7.24. 
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Table 7.2 Volume Fraction for Each Droplet Size Bin Specified in 
Example 4 


Bin Droplet Diameter 

Cumulative 

Volume 

Middle Right-Hand 

Volume 

Fraction 

Endpoint 

Fraction 



1 

4 

8 

0.0082 

0.0082 

2 

12 

16 

0.2119 

0.2037 

3 

20 

24 

0.7881 

0.5762 

4 

28 

32 

0.9918 

0.2037 

5 

36 

40 

1.0000 

0.0082 


When a droplet size distribution is input, the execution of the program 
is similar to when a single droplet size is specified. The difference is that 
the particle trajectory calculations must be performed for each droplet size 
specified. Since five discrete droplet sizes were specified in this example, 
the impingement limits and y 0 vs. s data were calculated for each of the 
droplet sizes. This will require significantly more computational time. 

This example will require user interaction when asked to preview the y 0 
vs. s data. This question will be asked for each droplet diameter, and the 
user will have the opportunity to remove points from each set of yo vs. s 
data. The procedure for removing the points was discussed in Example 1. 

Figures 7.25a-o show the plotted output created from the plot menu. 
Note that Figures 7.25b-f show the y 0 vs. s curves for each of the droplet 
diameters specified. These plots are made from a single selection of plot 
option 2. The user will be requested to input the axes limit for each set of 
y 0 vs. s data. The total local collection efficiency curve (Figure 7.25g) is 
calculated by su mming the contributions from each of the droplet diame- 
ters, as described in Section 4.2.5.I. The local collection efficiency for each 
droplet size is not stored by the program. 

Recall that in Example 1, an ice accretion was formed for this same icing 
condition except that a single droplet diameter of 20 microns was specified. 
A comparison of the local collection efficiency from the first timestep of 
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Examples 1 and 4 is shown in Figure 7.26. The local collection efficiency 
curves are similar near the stagnation point (s = 0.0). However, the region 
over which droplets impinge is much greater when a droplet distribution 
is specified because of the presence of the larger droplets. This result is 
similar to that concluded by previous investigators (Chang, Reference 13). 

The increased region of droplet impingement is also noticeable in com- 
parisons among the predicted ice accretion shapes. However, because of the 
small amount of mass present at the impingement limits, the effect on the 
total ice accretion shape is small. These results are not sufficient for deter- 
mining the total effect of using a multidispersed droplet size distribution 
to characterize an icing cloud as opposed to a monodispersed distribution. 
The effect of wider distributions with similar mass median droplet diam- 
eters must be evaluated. The effects must also be determined for various 
types of icing conditions and ice accretions. 

7.5 Example 5: Thermal Anti-icing 

The formulation of the thermodynamic equations in LEWICE allows the 
surface of the body to be a heat source or sink. However, since conduction 
through an existing ice layer is not modeled, the equations are valid only 
on the first timestep when ice is accreted on the clean airfoil. Also, since 
the melting of previously deposited ice from the surface of the body is 
not modeled, the equations are not applicable to a de-icing system. They 
are applicable, however, to the evaluation of a thermal anti-icing system. 
This example demonstrates the capability of LEWICE to determine the 
minimum thermal anti-icing heating requirements for maintaining an ice- 
free surface for a specified set of icing conditions. A discussion of the 
applicability of the thermodynamic model is presented in Appendix A. 

The NACA 0012 airfoil evaluated in this example has an electro-thermal 
anti-icing heater covering the first 20 percent of the airfoil, as shown in 
Figure 7.27. A uniform heat flux of 6000.0 W/m 2 was specified over the 
entire surface. Lewice is to be applied to determine whether this heat flux 
is sufficient to maintain an ice-free surface at the following icing condition: 

Angle of attack = 6.0 

Velocity = 128.41 m/sec 
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Static Temperature 
Static Pressure 
LWC 

Droplet Diameter 
Roughness Height 


= 265.0 K 
= 90748.0 Pa 
= 0.50 g/m 8 
= 20.0 microns 
= 0.0002 m 


The input file is shown in Figure 7.28. The heat flux from the surface (into 
the thermodynamic control volume) is input in namelist ICE through the 
variable QCOND. 

The icing time for this case was arbitrarily specified to be 120.0 sec so 
that any ice that did form would be easily visible on the plots. No user 
interaction, except the standard responses discussed in Examples 1 and 3, 
was required. As shown in Figure 7.29a, the heat input was not sufficient 
to maintain an ice-free surface. The ice accretion on the lower surface was 
formed as liquid water flowed off the heater and onto the colder airfoil 
surface. This result indicates that it is possible for ice to form aft of the 
heater but, since the shedding of water droplets from the surface are not 
modeled, the amount of ice may be over-predicted. The ice accretion on 
the upper surface was formed on the heater itself because the increased 
velocity (Figure 7.29d) over the upper surface caused a large decrease in 
the local static temperature (Figure 7.29e). It is also interesting to note 
the calculated surface temperature and convective heat transfer coefficient 
profiles over the surface (Figures 7.29g and h, respectively). The convective 
heat transfer was at a maximum on the upper surface which, when com- 
bined with the decreased local static temperature, contributed to the ice 
growth. The surface temperature profile shows that the temperature was 
greater than 0.0° C over a large portion of the heater, but was less than 0.0 
over a small portion of the upper surface. 

These results indicate that the heat flux should be increased in this 
region to maintain an ice-free surface. The actual minimum icing require- 
ment could be determined by incrementally increasing the heat input and 
re-calculating the ice accretion shape. Unfortunately, there was no exper- 
imental data to verify these calculations. This example was presented to 
demonstrate the capability of the thermal model. It is expected that this 
capability be further exercised to evaluate its accuracy. 
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NACA 0012 : 

EXAMPLE 1 


CS29Y 



ILIFT 3 1 



IPARA 3 1 



IFIRST* 3 



ISECND 3 3 



IPVOR 3 1 



INCLT* 0 



CLT 3 0.0 



ICHQRD 3 0 



CCL = 0.0 



IND 3 1 



ISOL= 0 



IPRINT* 0 



IFLLL 3 1 



LEND 



1 . 00C0300 

0 . 9899999 

0 . 9800000 

0.9000000 

0.8750000 

0.8500000 

0.7500000 

0 .7250000 

0 .7000000 

0 .6000000 

0 .5750000 

0.5500000 

0.9500000 

0 .9250000 

0 .9000000 

0 . 3000000 

0.2750000 

0.2500000 

0 .1500 000 

0.1250000 

0 .1000000 

0 . 0600000 

0 .0500000 

0 .0950000 

0 . 0250000 

0 . 0200000 

0 .0 150000 

0 . 0037500 

0 .0025000 

0 .0022500 

0 .00 12500 

0.0010000 

0 .0 008750 

0 . 0003750 

0 .0002500 

0 .0001250 

0 .0003750 

0 . 0005000 

0 . 0006250 

0 .00 12500 

0 .0015000 

0 .0017500 

0 .00 37500 

0 .0050000 

0 .0075000 

0 .0250000 

0 .0300000 

0 . 0350000 

0 . 0600000 

0 . 0700000 

-0 . 0800000 

0.1500000 

0 .1750000 

0 .2000000 

0 . 3 0 000 0 0 

0 . 3250000 

0 . 3500000 

0 .9500000 

0 .9750000 

0 . 5000000 

0 .6000000' 

0.6250000 

0 .6500000 

0.7500000 

0 .7750000 

0 .8000000 

0 . 90 00000 

0 . 9250000 

0 . 95000 0 0 

1 . 0000 000 

0.0000000 

0.0000000 

0 .0000000 

-0 .0018790 

-0.0036110 

-0 . 0 1 9 9080 

-0 .0179530 

-0 .0208880 

-0.0318550 

-0 .0399110 

-0 .036387 0 

-0 . 0959090 

-0 . 0979060 

-0 . 0997930 

-0 .0 56051 0 

-0 .0572930 

-0.0582620 

-0 . 060226 0 

-0 . 0600760 

-0.0596120 

-0 .0536 360 

-0.0507320 

-0.0970090 

-Q .0385350 

-0 .0356790 

-0.0390820 

-0 .0262300 

-0 .0237260 

-0 . 0207970 

- 0.0106990 

- 0 .0086200 

- 0.0081360 

-0.0058970 

-0.0051960 

-0 .0097660 

-0.0029950 

-0.0023560 

-0.0016320 

0.0029950 

0 .0039620 

0.0039310 

0 .0058970 

0 .0069850 

0.0070750 

0.0106990 

0.0123860 

0.0150780 

0.0262300 

0.0289620 

0.0309960 

0.0385350 

0.0910380 

0.0932520 

0.0536360 

0.0558790 

0 . 0575700 

0.0602260 

0.0600990 

0.0597070 

0.0560510 

0.0596980 

0.0531980 

0.0959090 

0.0937950 

0.0915850 

0.0318550 

0 . 0292210 

0.0265150 

0.0199080 

0.0117000 

0.0082510 

0 .0000000 

0 .0000000 

0.0000000 


ETRAJ1 

GEPS= 0 .9999999E-09 
DSHIFT* 0 . 2QE-02 
VEPS= 3.99999 9 9E- 0 3 
LCMB* 0 
Lcnp= o 

LEQfl 3 1 

LSi'n= o 

Li OR 3 1 
LXOR= 1 
N3DY = 1 
NEQ= 9 
NPL = 15 
NSEAR 3 50 
NSI = 1 

TIMSTP 3 0 . 9999999E-03 

LEND 

ET3AJ2 

CHORD- 0.30 

G= 0.0 

PIT 3 0.0 

PITDOT 3 0.0 

PRATK 3 0.0 

XORC= -9.0 

XSTOP* 0.50 

YOLIM 3 0 . 9999999E-09 

Y0P1AX* 0.50E-01 

Y0MIN* -0 . 50E-0 1 

YORC* 0 . 9999996E-01 

LEND 

CDIST 

FLUC* 1.0, 9*0.0 
DPD* 20.0, 9*0.0 
CFP* 1.0, 9*0.0 
£END 
DICE 

VINF* 129.960 
LWC= 0.50 
TAMB a 260.5998 
PAMB* 90798.0 
RH* 100.0 
DPHM- 20.0 
XKINIT 3 0.00035 
SEGTOL 3 1.50 


0 . 9700000 
0.8250000 
0.6750000 
0.5250000 
0.3750000 
0.2250000 
0 .0900000 
0 . 0900000 
0 .0100000 
0 .0020000 
0 .0007500 
0 .0000000 
0.0007500 
0.0020000 
0.0100000 
0 . 0900000 
0.0900000 
0.2250000 
0 .3750000 
0 . 5250000 
0.6750000 
0 .8250000 
0.9700000 
0.0000000 
-0 .0052390 
-0 . 0237390 
-0.0392810 
-0 . 0515600 
-0 .0590900 
-0.0587990 
-0.0952280 
-0.0323620 
-0.0172980 
- 0 . 0076230 
-0.0093630 
0 .0000000 
0.0093630 
0.0076230 
0.0172980 
0.0323620 
0.0952280 
0.0587990 
0.0590900 
0.0515600 
0.0392810 
0.0237390 
0.0052390 
0 . 0000000 


0 . 9500000 
0 .8000000 
0.6500000 
0 . 5000000 
0.3500000 
0 .2000000 
0 . 0800000 
0.0350000 
0 .0075000 
0 .0017500 
0 .0006250 
0.0001250 
0 .0008750 
0.0022500 
0 .0150000 
0 . 0950000 
0.1000000 
0.2500000 
0 .9000000 
0.5500000 
0 .7000000 
0.8500000 
0 . 9800000 
0 .0000000 
-0 .0082510 
-0 . 0265150 
-0 . 0915850 
-0.0531980 
-0 .0597070 
-0 .0575700 
-0.0932520 
-0.0309960 
-0 .0150780 
-0 . 007 0750 
-0.0039310 
0 .00 16320 
0.0097660 
0.0081360 
0.0207970 
0.0390820 
0.0970090 
0.0596120 
0.0582620 
0 .0997 930 
0 .0368870 
0 .0208880 
0.0036110 
0 .0000000 


0.9250000 
0 .7750000 
0 .6250000 
0 .9750000 
0.3250000 
0 . 1750000 
0 . 0700000 
0 . 0300000 
0 .0050000 
0 .0015000 
0 . 0005000 
0 . 0002500 
0.0010000 
0 .0025000 
0 . 0200000 
0 . 0500000 
0 . 1250000 
0 .2750000 
0 .9250000 
0.5750000 
0 . 7250000 
0 .8750000 
0 . 9899 999 
0 . QO000Q0 
-0 .0 1 1 7 0 0 0 
-0 . 0292210 
-0 . 0937950 
-0 . 0596 980 
-0 . 0600 990 
-0 .0558790 
-0 . 0910330 
-0.0289623 
-0 .0123860 
-O .0069850 
-0 . 0 0 39620 
0 .0023560 
0 .0051960 
0 .0086200 
0.0237260 
0.0356790 
0.0507320 
0 . 0600760 
0.0572930 
0 .097 9060 
0.0399110 
0 .01 7 9530 
0 .0018790 
0 . 0000000 


a. Input data file 


Figure 7.1: Input data for Example 1 


ORIGINAL PAGE IS 
OF POOR QUALITY 


6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 C 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
6 0 3 
1 1 3 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
6 0 9 
1 1 9 
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Time. » 0. 

ENTER DESIRED ICING TIME (SEC), (FIB. 01 

i?e. 

ENTER DESIRED TIME INCREMENT (SEC), (FI0.0) 

60 . 


are new plot options desired? iy/nj 
1/ 

available plot options 

0 - no plots 

1 - PARAMETER plots only 

2 - trajectory plots only 

3 - PARAMETER AND TRAJECTORY PLOTS 
ENTER PLOT OPTION (II) 

3 

THE POTENTIAL FLOW FIELD IS NOW BEING CALCULATED 


b. Interactive input 
Figure 7.1: Concluded 
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vElCCITY CHS S> 129.00 
TEMPERA TOPE <C> -12.60 
PRESSURE <<PA) 90.75 
"CM! CITY <T> 100.00 
lXC (3yy--3> 0.50 
CROP C r AK (MICRONS) 20. CO 
TIME (SEC) 0.00 


ORIGINAL PAGE IS 
OF POOR QUALITY 



Figure 7.2: y 0 vs s points calculated in subroutine COLLEC for Example 1 


THE CALCULATED SURFACE IMPACT DISTANCES (S> AND RELEASE POINTS ( YO ) 
ARE AS FOLLOWS , e 

; <0 S I Y0 S 


0.00706 

0.02715 

10 

- 0.00106 

- 0.00160 

O . 0 C 716 

e . 01944 

1 1 

- 0.00208 

- 0.00256 

0.00613 

0.01302 

12 

- 0.0031 1 

- 0.00385 

0.0091 1 

0.00064 

1 3 

- 0.00414 

- 0.00734 

0.02408 

0.00604 

14 

- 0 . 0051 7 

- 0.00983 

0.00309 

0.00370 

15 

- 0.00610 

- 0 . 01 376 

0 . 0 G 2 C 3 

0.00240 

16 

- 0.00722 

- 0.01886 

0.02100 

0.00154 

1 7 

- 0.00802 

- 0.03702 

- 0.00003 

- 0.00003 

18 

0.00000 

0.00000 


HOW -ANf TRAJECTORIES ARE TO BE DELETED BEFORE 
r hi. COt.LECTICN EFFICIENCY CALCULATION? < 1 2 J 


NACA 0012 : EXAMPLE I : TIME- 120.000 SEC 


STATIC TEMPERATURE (C) 

260. 

,55 

STATIC PRESSURE (PA) 

00748. 

,00 

VELOCITY (M/S) 

120. 

00 

LVC > G/M**3 ) 

0. 

50 

DROPLET DIAMETER (MICRONS) 

20. 

00 

:CRET ION DATA 

74 


STAGNATION POINT 


TRANSITION POINTS ( LOWER , UPPER ) 

72 


ICING LIMITS (LOWER. UPPER) 

43 


NUMBER OF POINTS 

147 


NUMBER OF SEGMENTS ADDED 

8 
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AVAILABLE PLOT OPTIONS 

0 - no plots 

1 - ICED GEOMETRY 

2 - Z VS S 

3 - BETA VS S 

4 - VE VS S 

5 - TE VS S 

6 - PE VS S 

7 - TSURF VS S 

8 - HTC VS S 

n _ y/v VR S 

10 - ICE DENSITY VS S 
I I - FFRAC VS S 
ENTER OPTION NUMBER (12) 


Figure 7.3: Interactive input and output for the second time step of Example 1. 
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S (1) 


a. Iced airfoil geometry d. Edge velocity vs. surface location 



b. Particle release position vs. surface 
impact distance 


o r 



S ( M ) 

e. Edge temperature vs. surface location 



S (Hi - ' f * > 

c. Local collection efficiency vs. surface f. Edge pressure vs. surface location 

location 


Figure 7.4: Icing parameter plots for the second timestep of Example 1. 
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(3) Si bO ~ " OlH 


ORIGINAL PAGE IS 
0F «X» QUAUry 


VELOCITY (M/S) 129.4-6 
TEMPERATURE (C) -12.60 
PRESSURE (KPA) 90.75 
HUMIDITY (*) 100.00 
LHC <G/Mhh3) 0.50 
DROP 01 AM (MICRON) 20.00 
TIME (SEC) 120.00 



-. 29 ~ r 7T\ TTcT - .' 1 o'o“t".c> c v "o'Tz: cTTI :.3 = 

s CM) 

Equilibrium surface temperature vs. 
surface location 


1 oco 

900 

BOO 

700 

2 t>00 
o 

W 500 
o 

- *00 
300 
200 
1 OC 


5 .75 TTTE r m -.'i* - .'3? -.do o.o? o.'i* 0.21 g .26 073 5 
s cm; 

j. Ice density vs. surface location 


BOO f 



S CM) 

h. Convective heat transfer coefficient 
vs. surface location 


. 0 0050p 
.00095 - 
. 000*0 - 
.00035- 
.00030- 
,00025- 
*. 00020 - 
.00015- 
.00010- 


. 0 D005p 

.ooogn 55 .. : g, 7TTV -.fry fl.fr ? C.m C. ' 2i C. : ?6 77 35 

S CM) 

i. Equivalent sand-grain roughness 
height vs. surface location 



t 

C.7 f- 


0 . o - 

.0 .'5 - • 

' 0 .- - 

0.3 - 

0.2 - 
I 

0. ! j- 


I 



0 *- c . - .07- - . 6 c e o TUI's 


D • £) 2 ** C • 0 ** 0 


k. Freezing fraction vs. surface location 


Figure 7.4: Concluded 
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ORIGINAL PAGE IS 
OF POOR QUALITY 

UMTRAMSrORHED COORDINATE DATA FOR BOOT ID - 1. NACA 0012 : EXAMPLE 1 


I 

xm 

VC I ) 

I 

XCI) 

Y( I) 

I 

XCI) 

Y( I ) 


1 .0 00000 0 

0 .0000000 

51 

0.0122522 

-0.0250171 

101 

0 .0339538 

0 .0331 726 

2 

0 . 9899997 

-0 . 0018790 

52 

0.0063893 

-0.0218130 

102 

0 .0391365 

0.0397692 

3 

0.9799998 

-0.0036110 

53 

0.0029776 

-0.0197319 

103 

0 . 0992957 

0.0362028 


0.9699997 

-0.0052390 

59 

-0.0006996 

-0.0170765 

109 

0 . 0999636 

0.0379372 

5 

0 . 99 99997 

-0.0082510 

55 

-0.0020565 

-0 .0 162655 

105 

0 . 0596575 

0 .0 397926 

6 

0 . 929 9998 

-0.0117000 

56 

-0.0039135 

-0 .0 159595 

106 

0 . 0698789 

0.0915216 

7 

0 .8999999 

-0.0199080 

57 

-0.0053179 

-’0.0127929 

107 

0.0799999 

0 . 09 32520 

8 

0 .879 9999 

-0 .0179530 

58 

-0.0053978 

-0 .0129796 

108 

0 . 0899999 

0 .0952280 

9 

0.8999998 

-0 . 0208880 

59 

-0.0053782 

-0.0119662 

109 

0.0999998 

0.0970090 

10 

0.8299999 

-0 . 0237390 

60 

-0.0059665 

-0 .0111969 

110 

0 . 1250000 

0.0507320 

11 

0 .7 999997 

-0.0265150 

61 

-0.0055827 

-0.0103023 

1 1 1 

0 . 1 9 999 9 9 

0 .0536360 

12 

0 .7 799997 

-0 . 0292210 

62 

-0 .0056817 

-0.0099199 

112 

0 . 1799999 

0 .0558790 

1 3 

0.7500000 

-0.0318550 

63 

-0.0057773 

-0 . 0059782 

113 

0 . 1999999 

0 .0575700 

1* 

0 . 729 9998 

-0 .0399110 

69 

-0.0058929 

-0 .0075095 

119 

0 .2299998 

0 . 0587990 

15 

0.6999997 

-0 .0368870 

65 

-0 .0059092 

-0.0069399 

115 

0 .25 0 0 0 0 0 

0.0596120 

1 -> 

0 .6 799998 

-0 . 0392810 

66 

-0 .0059506 

-0.0063720 

116 

0 . 2799997 

0 . 0600760 

17 

0.6999999 

-0 . 0915850 

67 

-0 .0060122 

-0 . 0057760 

117 

0.2999997 

0 . 0602260 

18 

0 . 6299999 

-0 . 0937950 

68 

-0 .0060896 

-0.0051370 

118 

0.3299999 

0 . 0600990 

19 

0 . 5999998 

-0 . 0959090 

69 

-0.0061991 

-0.0099373 

119 

0.3999998 

0 .0597070 

20 

0 . 5 79 9999 

-0 . 097 9060 

70 

-0.0063691 

-0 .0036290 

120 

0.3799999 

0.0590900 

21 

0 . 5999997 

-0 . 09 97930 

71 

-0.0069739 

-0 .0030971 

121 

0 . 3999999 

0 .0582620 

C L 

0 . 529 9997 

-0 .051 5600 

72 

-0.0065787 

-0.0029653 

122 

0.9299998 

0 . 0572930 

23 

0 . 5000000 

-0 .0531 980 

73 

-0.0066632 

-0 .0012327 

123 

0.9999997 

0 .0560510 

29 

0.9799995 

-0 .0 596980 

79 

-0.0067976 

0 .0000000 

129 

0.9799998 

0.0596980 

25 

0.9999997 

-0.0560510 

75 

-0.0066629 

0.0012326 

125 

0.5000000 

0 .0531980 

26 

0.9299998 

-0 . 0572930 

76 

-0.0065782 

0 .0029652 

126 

0 . 5299997 

0 .0515600 

27 

0 . 3999999 

-0.0582620 

77 

-0.0069731 

0.0030970 

127 

0 . 5999997 

0 .0997930 

28 

0 .379 9999 

-0.0590900 

78 

-0.0063680 

0.0036287 

128 

0 . 5799999 

0.0979060 

29 

0 . 3999998 

-0 . 05970 70 

79 

-0.0061976 

0.0099369 

129 

0 .5999998 

0 .0959090 

30 

0 . 3299999 

-0 . 0600 990 

80 

-0.0060878 

0 . 0051366 

130 

0 .6299999 

0 . 0937950 

31 

0 .2999997 

-0 . 0602260 

81 

-0.0060101 

0.0057755 

131 

0 .6999999 

0 . 0915850 

32 

0 .2799997 

-0 . 06 0 0760 

82 

-0.0059985 

0.0063713 

132 

0 .6799998 

0 . 0392810 

33 

0.2500000 

-0 .0596120 

83 

-0.0059070 

0.0069337 

133 

0.6999997 

0.0368870 

3 

0.2299998 

-0 . 0587990 

89 

-0.0058902 

0.0075087 

139 

0.7299998 

0 .0399110 

35 

0 . 1999999 

-0 . 0575700 

85 

-0.0057750 

0.0089779 

135 

0 .7500000 

0.0318550 

26 

0 .1799999 

-0.0558790 

86 

-0.0056799 

0.0099139 

136 

0.7799997 

0.0292210 

J 7 

3 . 1*99999 

-0 . 0536360 

87 

-0.0055800 

0.0103011 

137 

0 .7999997 

0.0265150 

38 

3 . 1 250 000 

-0 . 0507320 

88 

-0.0059638 

0.0111951 

138 

0.8299999 

0 . 0237 390 

2 9 

0 . 0999998 

-0 . 0970090 

89 

-0.0053769 

0 .0 1 1 9659 

139 

0 .8999998 

0 . 0208880 

90 

0 . 0899999 

-0.0952280 

90 

-0.0053963 

0.0129788 

190 

0.8799999 

0.0179530 

91 

0.0799999 

-0.0932520 

91 

-0.0053162 

0.0129922 

191 

0.8999999 

0.0199080 

92 

0 . 0698895 

-0.0919999 

92 

-0.0033869 

0.0159395 

192 

0 . 9299998 

0 .0 1 1 7000 

93 

0 .0596677 

-0 .0 397590 

93 

-0.0020270 

0.0162932 

193 

0 . 9999997 

0 .00825 1 0 

99 

0 . 0 9 99722 

-0 .0379086 

99 

-0.0006677 

0.0170519 

199 

0 . 9699997 

0.0052390 

95 

0 .099J028 

-0 . 036 181 3 

95 

0.0029863 

0.0197231 

195 

0 .9799998 

0.0036110 

96 

0.0391921 

-0.0397989 

96 

0 .0063991 

0.0218071 

196 

0 . 9899997 

0 .00 18790 

97 

0.0339577 

-0 .0331625 

97 

0.0122537 

0.0250150 

197 

1.0000000 

0 .0 0 0 0 0 0 0 

98 

0 . 028 7 098 

-0 .031976 7 

98 

0.0179189 

0.0275927 . 




99 

0.0233695 

-0.0296603 

99 

0.0233668 

0.0296667 




50 

0 .0179198 

-0 . 0275907 

100 

0.0287069 

0.0319896 




DOUGLAS A [PC PA 

FT COMPANY TWO- 

DIMENSIONAL 

POTENTIAL FLOW 

PROGRAM 





COMBINED 

FLOW 

NACA 0012 : EXAMPLE 1 






ALPHA = 

0.000000 

ALPHA O = 

0 .00 00 1 3 

HO. OF BODIES 1 




CL = 

0.000002 

CHORD = 

1 . 000000 

TOTAL 

ELEMENTS 196 




BCDY ID 

= 1 

NACA 0012 : EXAMPLE 1 

HO. OF 

ELEMENTS 196 




I 


i 

s 

VT 

CP 

J 

SIGMA 

VM 

1 

0.9950025 

-0.0009536 

0 . 0029715 

-0.8268182 

0 .316 3716 

1 

-0.013191 

0 . 0 0 0 009 

z 

0 . 985 0 022 

-0 .0 027 57 3 

0 .0079 085 

-0.8886936 

0.2103125 

2 

-0 . 013129 

0 . 0 0 0 009 

3 

0.97500 1-3 

-0 .009936 5 

0 .0123352 

-0 .9189090 

0 . 1556159 

3 

-0 . 0 126 7 0 

0 . 00 0 0 09 

9 

0 . 96 00095 

-0 . 0 067780 

0.0197099 

-0.9950629 

0. 1068562 

9 

-0 .0 1 2000 

0 .000009 

5 

0 . 9375099 

-0.0100112 

0.0307528 

-0.9680377 

0.0629031 

5 

-0.011207 

0 .0 0 0 0 09 

6 

0 . 9 I 2 5 C 2 5 

-0 . 0133289 

0 . 0930058 

-0 . 9839131 

0 .0319150 

6 

-0.010559 

0.000009 

7 

D .3875013 

-0.0169970 

0 . 0552962 

-0 . 9951219 

0 .0097339 

7 

-0.010190 

0 . 00000 3 

3 

0 . 8 6 25 C 1 1 

-0.0199323 

0 . 0679789 

-1.0093230 

-0 .0086691 

8 

-0 .00 990 7 

0 .000002 

9 

0.8375007 

-3.0223233 

0.0797051 

-1.0127668 

-0.0256953 

9 

-0.009760 

0.000001 

1 0 

0 . 8 1 25 0 0 5 

-0.0251359 

0 .0919275 

-1 . 0207739 

-0 . 091 9788 

10 

-0.009630 

0 . 0 0 0 0 0 1 

1 1 

0 .787 500 3 

-0 . 0278768 

0 . 1091959 

-1.0289805 

-0.0577717 

1 1 

-0.009516 

0 . 00000 1 

12 

0 . 762 50 08 

-0 .0305973 

0 . 1163605 

-1.0360298 

-0.0733576 

12 

-0.009389 

0 . 00000 1 

13 

0 . 7 3 7 500 7 

-0.0331927 

0 . 1285719 

-1.0933102 

-0 . 0889953 

13 

-0.009213 

0.000001 

19 

3.7125006 

-0.0356590 

0 . 1907789 

-1.0503035 

-0 . 1031370 

19 

-0 . 00 90 36 

0 . 0 0 0 0 0 1 

15 

0 .6375005 

-0 .0 380996 

0.1529815 

-1.0572052 

-0.1176825 

15 

-0.008853 

0 . 0 0 0 0 0 1 

16 

0 .66250 09 

-0 .0909999 

0 . 1651806 

-1.0690011 

-0 .1320982 

16 

-0.008626 

0 .0 00 0 0 1 

1 7 

0 6 57 50 08 

- 0 . 0927021 

0 . 1773756 

-1 . 0706520 

-0. 1962955 

17 

-0.008386 

0 . 0 0 0 0 0 1 

18 

0 . 6 125008 

-0.0998629 

0 . 1895665 

-1.0772190 

-0 . 1609009 

IS 

-0.008112 

0 . 00 0 0 0 1 

19 

0 . 58 7 50 0 7 

-0 . 0969188 

0.2017531 

-1.0836995 

-0.1793937 

19 

-0.007813 

0.000001 

20 

0 . 56250 06 

-0 . 0988691 

0.2139356 

-1.0900660 

-0 . 1882929 

20 

-0 . 00 797 7 

0 . 00000 1 

21 

0 .537500 5 

-0 . 0506920 

0.2261137 

-1.0963779 

-0.2020926 

21 

-0 .0 0 7127 

0.000001 

22 

0 .512500 9 

-0.0523956 

0.2382879 

-1.1027031 

-0.2159538 

22 

-0.006738 

0 .00000 1 

23 

0 . 98 7 500 9 

-0.0539657 

0.2509569 

-1.1089926 

-0.2298691 

23 

-0 .0 0631 1 

0 .000000 

29 

0 .96250 08 

-0.0553936 

0.2626229 

-1.1153297 

-0.2939609 

29 

-0.005851 

0 .00 0 0 0 0 

25 

0.9375007 

-0.0566678 

0.2797838 

-1.1216908 

-0.2581892 

25 

-0.005326 

0 . 0 0 0 0 0 0 

26 

0.9125007 

-0.0577752 

0.2869919 

-1.1280737 

-0.2725996 

26 

-0 . 009 755 

0.000000 

27 

0 . 3875006 

-0.0587010 

0.2990959 

-1.1395110 

-0.2871151 

27 

-0 .009 0 99 

0.000000 

28 

0 . 36250 05 

-0.0599260 

0.3112962 

-1.1908377 

-0.3015099 

28 

-0 . 00 335 1 

0 .000000 

29 

0 . 33750 03 

-0.0599308 

0.3233999 

-1.1970585 

-0.3157925 

29 

-0 .0 02529 

0 . 000000 

30 

0 .3 129 999 

-0.0601935 

0.3355910 

-1.1531639 

-0.3297863 

30 

-0.001583 

0.000000 

31 

0 . 2879 999 

-0 . 0601882 

0.3976869 

-1.1590567 

-0.3939119 

31 

-0.000526 

0 .0 0 0 0 00 

32 

0 .2629 990 

-0.0593856 

0.3598337 

-1.1697253 

-0.3565895 

32 

0 .0 0067 1 

-0 .000000 

33 

0.2379983 

-0 . 0592503 

0.3719836 

-1.1701090 

-0 . 3691925 

33 

0 . 002097 

-0.000000 

39 

0.2129971 

-0 . 0582363 

0.3891900 

-1.1799191 

-0.3809390 

39 

0 .003666 

0.000000 

35 

0 . 18 79 956 

-0 . 0567876 

0.3963069 

-1.1788921 

-0 . 3897858 

35 

0 .0 05558 

-0 . 000002 

36 

0 . 1629 931 

-0 . 0598323 

0.9089910 

-1 . 1819556 

-0 .3958368 

36 

0 .007853 

-0 .000002 

3* 

0 . 1 379893 

-0.0522799 

0.9207022 

-1 . 1815958 

-0.3961678 

37 

0 .010653 

-0 .000009 

38 

0 . 1 1 298 3 1 

-0 . 0989800 

0.9329569 

-1.1766272 

-0 . 3899509 

38 

0.019293 

-0 .000009 

39 

0.0999959 

-0 . 0961379 

0.9915639 

-1.1659998 

-0 . 3599389 

39 

0 .017908 

-0 . 0 0 0 009 

90 

0 . 089 999 9 

-0 . 0992900 

0.9965072 

-1.1386397 

-0.2969888 

90 

0 . 019990 

-0.000003 

91 

0.0799938 

-0.0923663 

0.9519770 

-1.1286559 

-0.2738628 

91 

0 .019238 

-0 .000003 


Figure 7.5: Printed output for the second time step of Example 1 
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DOUGLAS AIRCRAFT COMPANY TWO-DIMENSIONAL POTENTIAL FLOW PROGRAM 


COMBINED 
ALTHA = 

CL = 
BODY ID = 


FLOW MAC A 0012 : EXAMPLE 1 

0.000000 ALPHA O = 0.000013 

0 . 000002 CHORD = 1 .000000 

1 NACA 0012 : EXAMPLE 1 


NO. OF BODIES 1 
TOTAL ELEMENTS 146 
NO. OF ELEMENTS 146 


I 

X 

Y 

A 2 

0 . 064 7 76 1 

- 0 . 04 0626 7 

43 

3 . 3545632 

-0 .0 386 1 06 

uc. 

0 . 04 6 8354 

-0 .0 366 035 

<* 5 

0 . 04 17 165 

-0 .03548 50 

46 

0 . 0 36 54 6 1 

-0.0339678 

47 

0 .0313305 

-0.0323297 

48 

0.0260337 

-0 .0 30 5865 

4? 

0 .0206328 

-0 . 0286 31 3 

50 

0 .0 1 50664 

- 0 . 026 3229 

51 

0 .0 0 92 924 

-0.0234670 

52 

0 . 0046 7 02 

-0 . 0207942 

53 

0.0011390 

-0 . 0 184042 

54 

-0 . 00 1 37 7 9 

-0.0166708 

« ^ 

-C .0^27 350 

-0 .0 1 58600 

56 

-0 . 00458 92 

-0.0143968 

57 

-0.0053327 

-0.0127363 

55 

-0 .005 36 30 

-0 . 0 1 22229 

59 

-0 . 0 054 18 1 

-0 .0 1 1 5558 

60 

- 0 .0 0 55246 

-0 . 0 1 07243 

6 1 

-0 . 0 056 339 

-0.0098588 

62 

-0.0057255 

-0 . 0 089465 

6 3 

- 0 .0 358 348 

-0 . 007 9938 

64 

-0.0059008 

-0 .00 7221 9 

65 

-C . 0059274 

-0 . 0066530 

66 

- 0 . 9 0 5 9 7 9 7 

-0 . 00607 38 

67 

-0 . 004 048 9 

-0 . 00 5456 3 

62 

-0.0061407 

-0 .0047866 

69 

- 0 .006284 1 

-0 . 0 040 33 1 

: o 

-0 .0 06421 5 

-0.0033380 

7 1 

-0 . 006 526 3 

-0 .0 027562 

72 

-0.0066215 

- 0 . 00 18490 

73 

~C . 0 06 7 0 5 9 

-Q . 0006 164 

74 

-0.0067052 

0 . 0 0 06 16 3 

75 

-0 . 0 066 20 5 

0 . 0 0 18489 

76 

-0 .0 06 5258 

0 .0027 56 1 

77 

-0.0064006 

0 . 0 03337 9 

72 

- 0 . 0 06 2828 

0 . 0040328 

79 

-0 .006 1 390 

0 .0047862 

80 

-0 .0 06 047 0 

0.0054558 

81 

-0 .0059776 

0 .0060732 

82 

-0 .0 0 59253 

0 .0066523 

83 

-0 . 0058486 

0.0072212 

84 

-0 .0058326 

0 . 007 9930 

85 

-0 .0 057272 

0 . 0089456 

86 

-0.0056314 

0 .0098577 

87 

-0 . 0055219 

0.0107231 

88 

-0 .00 54 158 

0 .01 15548 

89 

-0.0053613 

0.0122221 

90 

-0.0053312 

0.0127355 

91 

-0 .0 045729 

0 .0 143884 

92 

-0 . 0027 06 7 

0 .0 1 58 388 

93 

-0.0013472 

0.0166473 

94 

0.0011593 

0 .0 18 387 5 

95 

0.0046765 

0 . 020 7875 

96 

0 . 0 092954 

0 . 02346 3 1 

97 

3 . 0 1 5 066 7 

0.0263228 

98 

0 . 0206 308 

0 .0286356 

99 

0 . 026 0 304 

0 .0 305937 

100 

3 . 0313269 

0 .0323384 

101 

0 .0365414 

0.0339804 

102 

0.0417105 

0.0355036 

103 

0 . 0468 7 75 

0 .0368288 

1 C 4 

0 .0 5455 34 

0 .0 386457 

105 

0.0647681 

0.0406571 

106 

0 . 074 9406 

0.0423782 

107 

0 . 084 9999 

0 .0442400 

108 

0 .0 94 9959 

0.0461379 

109 

0 .1 12483 1 

0 .0489800 

110 

0.1374893 

0.0522749 

111 

0 .1 62493 1 

0.0548323 

112 

0 . 18 74 956 

0.0567876 

113 

0.2124971 

0.0582363 

114 

0 .2374 983 

0.0592503 

i 1 5 

0.2624990 

0.0598856 

116 

0.2874994 

0.0601882 

117 

0 . 3124999 

0 . 0601935 

118 

0 .337 5003 

0.0599308 

119 

0 . 3625005 

0.0594260 

120 

0 . 387 5006 

0 .0587 0 10 

121 

0 .4 125007 

0 .0577752 

122 

0 .4 375007 

0 .0566678 

123 

0.4625008 

0 .0553936 


S 

VT 

CP 

0.4564384 

-1 . 1665764 

-0 . 36 08999 

0.4615475 

-1.1738958 

-0.3750308 

0.4653795 

-1 .1 783762 

-0 . 3885698 

0.4679711 

-1 .1814480 

-0 . 3958187 

0.4705892 

-1 .16 33625 

-0 . 353*117 

0.4732450 

-1 . 149594 3 

-0.3215666 

0.4759541 

-1.1521044 

-0.3273439 

0.4787447 

- 1 . 14 99968 

-0 . 3224 91 6 

0.4816726 

-1 . 140 1 50 1 

-0 .2999420 

0 . 4848031 

-1.0972996 

-0 . 204066 3 

0.4873973 

-1.0642109 

-0.1325445 

0 .4894699 

-0 . 9627361 

0 .07 31 393 

0.4909557 

-0.9737462 

0 .0518184 

0.4917236 

-1.1018343 

-0.2140379 

0.4928801 

-1.2881260 

-0.6592684 

0.4937775 

-0.8853139 

0.2162194 

0.4940273 

-0 .7 1 047 51 

0.4952251 

0.4943524 

-0.6279104 

0.6057286 

0.4947596 

-0.5720486 

0.6727604 

0.4951834 

-0 .5080 939 

0.7418406 

0.4956290 

-0.4414389 

0.8051318 

0.4960946 

-0.4165254 

0.8265067 

0.4964712 

-0.3197475 

0.8977616 

0.4967479 

-0 . 3052155 

0 . 9068435 

0.4970304 

-0.2758528 

0.9239053 

0.4973322 

-0.2470832 

0 . 9389499 

0 .4976606 

-0.2164860 

0.9531339 

0.4980332 

-0.2118799 

0 .9551 06 9 

0.4983774 

-0 .1746964 

0.9694812 

0.4986646 

-0 . 17836 13 

0. 9651873 

0 .4991083 

-0 .093594 3 

0 . 9912401 

0.4997085 

-0 . 06 07 150 

0 .996 31 38 

0.5003088 

0 .06 02920 

0 .9963649 

0.5009090 

0 .09447 92 

0.9910737 

0.5013527 

0.1792508 

0 . 9678692 

0 .5016398 

0 .1 75031 1 

0 . 96 93642 

0 .5019841 

0 .2121464 

0 . 9549940 

0.5023567 

0.2167311 

0.9530277 

0.5026851 

0.2474253 

0.9387807 

0 . 5029869 

0.2760307 

0.9235071 

0.5032694 

0 . 3056121 

0.9066013 

0 . 5035461 

0 . 3200547 

0.8975651 

0.5039227 

0.4169908 

0.8261187 

0.5043883 

0.4419479 

0.8046821 

0.5048338 

0.5088874 

0.7410337 

0.5052575 

0.5725474 

0.6721895 

0.5056649 

0.6287807 

0.6046349 

0.5059901 

0.7130577 

0.4915488 

0.5062399 

0.8891031 

0.2094955 

0.5071375 

1.2907629 

-0.6660681 

0.5082944 

1.0958996 

-0.2009954 

0.5090628 

0.9688777 

0.0612761 

0.5105464 

0.9584563 

0.0813615 

0.5126163 

1.0659027 

-0.1361485 

0.5152097 

1.0972567 

-0.2039719 

0.5183399 

1.1400375 

-0.2996855 

0.5212675 

1.1501760 

-0.3229046 

0.5240582 

1.1523075 

-0.3278122 

0.5267673 

1 . 1492682 

-0.3208170 

0.5294232 

1 . 1631508 

-0.3529196 

0.5320413 

1.1817484 

-0.3965292 

0.5346329 

1.1791773 

-0.3904591 

0.5384650 

1.1760740 

-0.3831501 

0.5435247 

1.1687498 

-0.3659754 

0.5485370 

1.1282282 

-0.2728957 

0.5535073 

1.1373749 

-0.2936211 

0.5584505 

1 . 1656160 

-0.3586607 

0.5670580 

1.1764851 

-0.3841162 

0.5793122 

1.1815319 

-0.3960171 

0.5915234 

1.1814203 

-0.3957539 

0.6037076 

1.1785692 

-0.3897324 

0.6158745 

1 , 17490 39 

-0 . 3803988 

0.6280306 

1.1700935 

-0 . 3691 187 

0.6401806 

1.1647167 

-0.3565645 

0.6523275 

1.1590490 

-0.3433943 

0.6644735 

1.1531582 

-0.3297729 

0.6766200 

1 . 1470537 

-0.3157320 

0.6887682 

1.1408339 

-0.3015013 

0.7009190 

1.1345081 

-0.2871084 

0.7130730 

1.1250708 

-0.2725430 

0.7252306 

1 .1216888 

-0.2581854 

0.7373920 

1.1153278 

-0.2439556 


ORIGINAL PAGE IS 
OF POOR QUALITY 


J 

SIGMA 

VM 

42 

0.011522 

-0.000002 

43 

0.019424 

-0.000002 

44 

0 . 0 18 50 3 

-0.000002 

45 

0 .0235 5 9 

- 0 .0 0 0 0 0 1 

46 

0 .026 6 92 

-0.000002 

47 

0 . 027 1 36 

-0.000002 

48 

0 .0272 43 

-0 . 0 0 0 00 1 

49 

0 . 0 3 1 244 

-0.000001 

50 

0 .0 34 985 

-0.000001 

51 

0 .042964 

-0.000000 

52 

0 . 04 3626 

-0.000000 

53 

0 . 0490 18 

0.000000 

54 

0 . 02 3584 

0 .000 002 

55 

0 .009558 

0 .000003 

56 

0 . 043306 

0 .000004 

57 

0 . 126 043 

-0 . 000032 

58 

0 . 128969 

-0.000030 

59 

0 . 128295 

- 0 . 0 0 0026 

60 

0 . 123 520 

-0.000018 

61 

0 . 1248 37 

-0.000024 

62 

0 . 12766 1 

-0 .0 0 0 025 

63 

0 .1256 30 

-0.000024 

64 

0 .131645 

-0.000032 

65 

0 . 13647 7 

-0.000023 

66 

0.135210 

-0.000023 

67 

0.134670 

-0 .00 0 022 

68 

0 . 1 33490 

-0.000015 

69 

0 .1 24528 

-0.000002 

70 

0 .11 9806 

-0.000010 

71 

0 . 1 1 4872 

-0.000008 

72 

0.118039 

-0.000025 

73 

0 . 1 16 084 

-0.000027 

74 

0 . 1 1620 9 

-0 . 00 0 026 

75 

0 .118143 

-0.000027 

76 

0 .1 147 37 

-0.000007 

77 

0 .1 1 97 28 

-0.000008 

78 

0 . 1 244 92 

- 0.000003 

79 

0 .1 33458 

-0 . 0 0 0 0 1 5 

80 

0.134660 

-0 .0 0002 1 

81 

0.135213 

-0 .000023 

82 

0 . 136443 

-0 .000026 

83 

0 .131615 

-0 .0 00026 

84 

0 . 12557 5 

-0.000019 

85 

0 . 127589 

-0 .000023 

86 

0.124722 

-0 . 00002 1 

87 

0.123468 

-0.000018 

88 

0 .128303 

-0 .000026 

89 

0 .128798 

-0 .000025 

90 

0.125787 

-0.000024 

91 

0 . 04 1 929 

0 . 00 0004 

92 

0 .009722 

0 .00000 1 

93 

0.023746 

-0.000001 

94 

0.050142 

-0 .0 0000 1 

95 

0.043791 

-0 . 0000 02 

96 

0 . 043104 

-0.000002 

97 

0.035132 

-0 .000002 

98 

0.031395 

-0 . 00000 3 

99 

0 . 027300 

-0.000003 

100 

0.027225 

-0.000003 

101 

0.026901 

-0 . 000003 

102 

0.023755 

-0.000003 

103 

0.018736 

-0 .0 0 0004 

104 

0.019568 

-0 . 000004 

105 

0.011189 

-0.000003 

106 

0.013823 

-0 . 0000 04 

107 

0.019942 

-0 .000 0 0 3 

108 

0 .0 1 7405 

-0.000004 

109 

0.014238 

-0 . 0000 0 3 

110 

0 .0 1 0649 

-0 . 000 002 

111 

0 .00 7848 

-0.000002 

112 

0.005554 

-0 .0 000 02 

113 

0.003662 

0.000000 

114 

0 .0 02044 

0 . 000 00 0 

115 

0 .000669 

0 . 0000 00 

116 

-0 .000528 

0 . 00 000 0 

117 

-0.001585 

0.000000 

118 

-0.002526 

0.000000 

119 

-0 .003353 

0.000000 

120 

-0.004101 

0.000000 

121 

-0.004757 

0 .000000 

122 

-0 .005328 

0.000000 

123 

-0.005852 

0 .00000 0 
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DOUGLAS 

A TPCRAFT COMPANY 

TWO-DIMENSIONAL POTENTIAL 

FLOW 

PROGRAM 


COMBINED 

FLOW MAC A 

0012 

; EXAMPLE 

1 



ORIGINAL PAGE IS 

ALPHA = 

0 .003000 

ALPHA O = 

0 .00 00 1 3 

HO. OF BODIES 1 

OF POOR QUALITY 

CL = 

0 . 000002 


CHORD = 

1 .000000 

TOTAL 

ELEMENTS 146 


BODY ID = 

1 NACA 

0012 

: EXAMPLE 

1 

MO. OF 

ELEMENTS 146 



I 

X 

Y 

S 

VT 

CP 

J 

SIGMA 

VN 

124 

0 . 48 7500 9 

0.0539657 

0.7495575 

1 .1089916 

-0.2298622 

124 

-0.006312 

0 .0 0 00 0 0 

125 

0.5125009 

0.0523956 

0 .7617270 

1 .1027021 

-0.2159519 

125 

-0 . 0 0 6 7 3 9 

0 . 000001 

126 

0 .5375005 

0 .0506920 

0 .7739007 

1.0963764 

-0.2020407 

126 

- 0 . 00 7 1 28 

0 .000001 

127 

0 . 5625006 

0 . 0488641 

0 . 7860789 

1.0900640 

-0.1882391 

127 

-0 . 00 74 78 

0 . 000001 

128 

0 . 5875007 

0.0469188 

0.7982613 

1 . 0836935 

-0.1743908 

123 

-0.007814 

0 .000001 

129 

0.6125008 

0 . 0448624 

0.8104479 

1.0772171 

-0. 1603966 

129 

-0 . 0 08 1 1 2 

0 . 000001 

130 

0 .6375008 

0 . 0427021 

0.8226388 

1.0706501 

-0 . 1462908 

130 

-0 . 008386 

0 .000001 

131 

0 .6625004 

0 . 0404444 

0.8348338 

1.0640001 

-0.1320953 

131 

-0 . 0086 27 

0.000001 

132 

0 . 6875005 

0 .0380946 

0.8470329 

1.0572042 

-0 . 1 1 76805 

132 

-0 . 0 08853 

0 . 000001 

133 

0.7125006 

0 . 0356590 

0.8592360 

1 .0503025 

-0.1031351 

133 

-0 . 0 0 90 36 

0 .000001 

134 

0 .7375007 

0 . 0331427 

0.8714430 

1.0433102 

-0 . 0884953 

134 

-0.009214 

0 .000001 

135 

0 .7625008 

0 .0305473 

0.8836539 

1.0360289 

-0 .0733557 

135 

-0 .009384 

0 . OOOCOl 

136 

0.7875003 

0 . 0278768 

0.8958685 

1 . 0284805 

-0 . 057771 7 

136 

-0.009517 

0 .0 0 0 0 0 1 

137 

0 .8125005 

0 .0251359 

0 . 9080870 

1 . 0207729 

-0 . 0419769 

137 

-0 . 0 0 96 3 1 

0 .0 0 0 0 0 1 

138 

0 .8375007 

0 . 0223233 

0 . 9203093 

1.0127659 

-0 .0256 939 

138 

-0.009760 

0 . 00000 1 

139 

0 .86250 1 1 

0 .0194323 

0 . 9325360 

1 .0043221 

-0 .0086622 

139 

-0 .0 09907 

0 . 00000 1 

140 

0 .88750 18 

0.0164470 

0 . 9447682 

0.9951208 

0 .0 0 97346 

140 

-0.010141 

0 .000001 

141 

0 . 9125025 

0 . 0133284 

0 . 9570085 

0.9839124 

0 .0319164 

141 

-0.010560 

0 .00000 1 

142 

0 . 9375044 

0.0100112 

0.9692615 

0 . 9680369 

0.0629046 

142 

-0 .011207 

0 .000004 

143 

3.9600045 

0 . 0067 780 

0 . 9803048 

0.9450626 

0 .1 068568 

143 

-0.012000 

0 .000005 

144 

0 . 9750013 

0 . 0044365 

0 . 9876790 

0.9189037 

0 .1556160 

144 

-0.012670 

0 . 000004 

145 

0 . 9850022 

0 .0027573 

0 . 9926056 

0.8886423 

0.2103150 

145 

-0.013129 

0 . 00000 5 

146 

0 . 9950 025 

0 .0009536 

0 . 9975426 

0.8268158 

0.3163756 

146 

-0 .013141 

0 .000004 


INTEGRATED VALUES 


CY = 

0 .0 0003 

CX = 

-0 .00106 



CL - 

0.00 0 0 3 

CD = 

-0.00106 

CM = 

0 . 00000 


PARABOLIC INTEGRATION 


INTEGRATED VALUES 

CY - 0.00002 CX = -0.00030 

CL - 0.00002 CD = -0.00030 CM = 0.00000 

TOTAL CM = 0.00000 

TOTAL Ctl = 0.0 000 0 (PARABOLIC) 


THE PARTICLES ARE RELEASED FROM X = -1.26000E 00 
WHICH IS OBTAINED AT THE 1 LOOP OF 50 LOOPS 

GEOMETRY CHARACTERISTICS : 

LEADING EDGE (X,Y) 

TRAILING EDGE (X,Y) 

THICKNESS 
CHORD 

ANGLE CF ATTACK 

UPPER BCUNDAPY 
LOWER BOUNDARY 

PARTICLE TRAJECTORY DATA: 


THE PARTICLES ARE RELEASED IN EQUILIBRIUM WITH THE AIR 


PARTICLE 

INITIAL 

INITIAL PARTICLE PITCH PIT 

GRAVT 

ERROR 


iH a:\cTER 

vx 

VY 

AOA 

ANGLE 

DOT 

CONST 

CRITERIA 


(MICRONS) 

(M/S) 

(M/S) (DEGREES) (DEGREES) (DEG/SEC) (M/S**2) 


20.00 

0 . 00 

0.00 0.00 

0.00 

1 0.00 

0 .00 

5. 00E-05 


THE PARTICLES OF SIZE 

2 0 . 0 OCOMTAIN 1.0000 OF THE TOTAL MASS 




xo 

Y0 



XP 

YP 

S 

DT 

NSTP 

-1.2599993 

0 .0090961 

OUT 

OF RANGE 

0.0831582 

0.0226937 


4.6785E-05 

101 

-1.2599993 

-0 .0090835 

OUT 

OF RANGE 

0.0828083 

-0.0226548 


5.2456E-05 

97 

Y0MAX* 9.09olE-03 Y0MIN= 

-9.0835E-03 






-1.2599993 

0 .00 00063 

HIT 

BODY AT 

-0.0040239 

0.0000050 

0.0000050 

8 . 12Q6E-06 

41 

-1.2599993 

0.0045512 

HIT 

BODY AT 

-0.0009029 

0.0054289 

0.0082756 

2.8464E-06 

68 

-1.2599993 

0 .0068236 

HIT 

BODY AT 

0.0061579 

0.0088386 

0.0162043 

9.2295E-06 

65 

-1.2599993 

0 .0079599 

HIT 

BODY AT 

0.0166598 

0.0117465 

0.0271505 

6.7323E-06 

80 

-1.2599993 

0 .0085280 

OUT 

OF RANGE 

0.0813009 

0.0221908 


4 . 7782E-0 5 

95 

-1.2599993 

0 .0 082439 

OUT 

OF RANGE 

0.0829488 

0.0222186 


4.7512E-05 

111 

-1.2599993 

0 . 0081019 

OUT 

OF RANGE 

0.0871353 

0.0226163 


6 . 1404E-05 

116 

-1.2599993 

0 .0080309 

OUT 

OF RANGE 

0.0854309 

0.0225082 


4.4515E-05 

119 

-1.2599993 

-0 . 0005618 

HIT 

BODY AT 

-0.0039709 

-0.0006448 

-0.0007154 

6.9591E-06 

45 

-1.2599993 

-0 .0048227 

HIT 

BODY AT 

-0.0003887 

-0.0057987 

-0.0089134 

3. 1107E-06 

56 

-1.2599993 

-0 . 0069531 

HIT 

BODY AT 

0.0069689 

-0 .0091094 

-0.0170804 

1 . 1028E-05 

58 

-1 .2599993 

-0 . 0080183 

HIT 

BODY AT 

0.0263552 

-0.0135138 

-0.0370157 

9. 9582E-06 

91 

-1.2599993 

-0.0085509 

OUT 

OF RANGE 

0.0868925 

-0.0228065 


4.9236E-05 

87 

-1.2599993 

-0.0082846 

OUT 

OF RANGE 

0.0823606 

-0.0221761 


4 . 0522E-05 

108 

-1 .2599995 

-0.0081514 

OUT 

OF RANGE 

0 . 0851743 

-0.0224173 


5. 1358E-05 

111 

-1 .2599993 

-0 . 0080849 

OUT 

OF RANGE 

0 . 0836582 

-0.0222142 


4 .8017E-05 

92 

UPPER SURFACE LIMIT 


LOWER SURFACE LIMIT 





YOU 

SU 


Y0L 

SL 





0.7960E-02 

0 .2715E 

-01 

-0 .8018E-02 

-0.3702E- 

01 




-l .2599993 

0.0071610 

HIT 

BODY AT 

0 .0082593 

0.0095326 

0 .0184434 

9. 592 9E-06 

71 

-1.2599993 

0 .0061338 

HIT 

BODY AT 

0.0032192 

0.0076855 

0.0130234 

7 .8691E-06 

68 





Figure 7.5: continued 




-4 . Q243E-0 3 2 . 07 15E-07 

3.0200E-01 0.0000 

4 . 0136E-02 
3 . 0000E-01 
0. 0000 

2 . 2075E-02 
-2 .2075E-02 
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-1 .2599993 0.0051066 HIT BODY AT 0.0002541 
-1.25 099 9 3 0. 0040775 HIT BODY AT -0.0014151 
-1.2599093 0.0030523 HIT BODY AT -0.0036111 
-1.2579975 0.0020251 HIT BODY AT -0.0037444 
-1.2599993 0.0009940 HIT BODY AT -0.0034594 
-1.2599993 -0.0000292 HIT BODY AT -0.0040216 
-1.2599993 -0.0010564 HIT BODY AT -0.0034440 
-1.2599993 -0.0020335 HIT BODY AT -0.0037397 
-1.2599993 -0.0031107 HIT BODY AT -0.0036065 
-1.2599993 -0.0041379 HIT BODY AT -0.0016157 
-1.2599993 -0.0051650 HIT BODY AT 0.0004149 
-1.2597993 -0.0061722 HIT BODY AT 0.0034358 
-1.2599793 -0.0072194 HIT BODY AT 0.0046541 


0.0061759 0.0096382 4.0084E-06 54 

0.0044164 0.0069433 2.4937E-06 52 

0.0035133 0.0037547 6.9544E-Q6 52 

0.0023366 0.0024918 1.0324E-05 56 

0.0011532 0.0015363 3.0499E-06 50 

-0.0000330 -0.0000332 8.1260E-06 39 

-0.0012178 -0.0016015 8.0285E-06 46 

-0.0024063 -0.0025623 6.3720E-06 56 

-0.0035746 -0.0038509 4.9189E-06 61 

-0.0044994 -0.0073419 2.7035E-06 59 

-0.0062647 -0.0098281 4.3314E-06 60 

-0.0077756 -0.0132626 9.5905E-Q6 63 

-0.0096553 -0.0188647 7.5941E-06 55 


YO VS 5 DAT.’. FOR DROPLET DIAMETER" 20.0 0000 MICRONS 
1' POINTS HAVE BEEN CALCULATED 



S 

YO 

i 

0.027150 

0.007960 

2 

0 . 018443 

0 . 0 0 7 1 6 1 

3 

0 . 0 1 3C23 

0 . C06 1 34 


P . 9G9638 

0 0 0 5 1 0 7 

5 

0 . G 0 6 9 4 3 

0 . 004 0 7 9 

6 

7.003785 

0 003052 

7 

0 . 002492 

0 . 002 025 

8 

0001536 

0 . 000998 

9 

-0 . 000033 

-0 . 000029 

1 0 

- 0 0 0 1 6 0 2 

-0 . 001056 

1 1 

0.002562 

-0 . 002044 

1 2 

-0 .003451 

-0.003111 

1 3 

-0.007342 

-0.004134 

14 

- 0 . 00 7428 

-0 . 005165 

15 

-0.013263 

-0.006192 

16 

-0 . 0 18 565 

-0 .0072 1 9 

17 

-0.037016 

-0 . 004014 


CALCULATED LOCAL 

COLLECTION EFFICIENCY 

FOR DROPLET 

DIAMETER* 20. 

.00000 MICRONS 

SEG 

5 

BETA 

SEG 

S 

BETA 

1 

-0 . 311805 

0 . 0 0000 0 

51 

-0 .011512 

0 .3 12600 

2 

-0.307626 

0 . 0 0000 0 

52 

-0 . 009732 

0 .334206 

3 

-9.304559 

0 . 000 0 00 

53 

-0.008418 

0 . 332866 

4 

-0.299981 

0 .00 0900 

54 

-0 .007 550 

0 . 3320 1 5 

5 

-0.293141 

0 . ocoooo 

55 

-0 .0068 94 

0.373618 

6 

-0 . 28 5557 

0.000000 

56 

-0 . 0057 36 

0 .4 086 03 


-0.277985 

0.000000 

57 

-0.004476 

0.446685 

8 

-0 . 2 704 22 

0.000000 

58 

-0 .00 3842 

0 .6721 9 7 

9 

-0 .26256 5 

0.000000 

59 

-0 . C03/04 

0 . 7020 0 l 

1 0 

-9.255312 

9.000000 

60 

-0.003433 

0 . 7 4 8 9 C 6 

; i 

- 0 2 4 7 7 6 1 

0.000000 

61 

-0.003210 

0 .80 9049 


-0.240213 

0.000000 

62 

-0.002921 

0.871553 

1 ; 

-9 .23266 7 

0.000000 

63 

-0.002613 

0 . 9 38362 

1 4 

- 0 . 225 1 22 

0 . 0 0 0 0 00 

64 

-0 .0 0232 1 

0 .850 024 

1 5 

-0.217580 

0.000000 

65 

-0 . 0G2153 

0 .8 37777 

16 

-0.210040 

0 . 000 000 

66 

-0 .00204 0 

0.829479 

1 7 

- 9 .202502 

0.000000 

67 

-0.001904 

0.819519 

iZ 

-0.194966 

0.000000 

63 

-0.001772 

0 .809 902 

1 9 

- 0 . 1 8 7 4 3 3 

0.000000 

69 

-0 .0 0 1 5 98 

0 .857 995 

20 

-9.179901 

0.000000 

70 

-0 . 0 0 1 38 1 

0.833184 

2 1 

-0.172371 

0 . 0 0 0 0 0 0 

71 

-0.001112 

0.802574 

22 

-0 . 164844 

0 . 0 0 0 00 0 

72 

-0 . 000728 

0 .758828 

23 

-0 . 1 57 3 14 

0 . 0 0000 0 

73 

-0 .000254 

0.704716 

24 

-0 . 149794 

0 .0 00 000 

74 

0 . 00 0 254 

0 .712874 

25 

-0 . 1 42272 

0 . oooooo 

75 

0 .000 729 

0.767481 

26 

-9 . 1 34750 

0.000000 

76 

0.001113 

0.811653 

27 

-0 . 127 230 

0.000000 

77 

0.001331 

0 . 842557 

24 

-0.119710 

0 .00 000 0 

78 

0 .00 1 599 

0.802467 

29 

-0.112190 

0 .000 00 0 

79 

0 . 00 1 773 

0 .815026 

30 

-0 . 104669 

0 .0 00000 

80 

0 .0 0 1 904 

0 .824530 

3 1 

-0 . 0 97 146 

0 .000000 

81 

0 .002040 

0.834363 

32 

-0 . 049620 

0 .000000 

82 

0.002153 

0.842550 

33 

-0 . 082059 

0 .00000 0 

83 

0.002321 

0 . 854658 

34 

-0 .074550 

0.000000 

84 

0.002612 

0 . 924 965 

35 

-0 . 067 000 

0 .000000 

85 

0 .002921 

0.857384 

36 

-0 . 0 59432 

0 . oooooo 

86 

0.003210 

0.794160 

37 

-0.051335 

0.000000 

57 

0 .003457 

0.733540 

34 

-0 . 0442 18 

o . oooooo 

88 

0.003702 

0.686370 

39 

-0 .038849 

o . oooooo 

89 

0.003840 

0.500825 

4 0 

-0 . 0 3577 0 

0 .oooooo 

90 

0.004483 

0 .473752 

4 1 

-0.032716 

0 .oooooo 

91 

0.005751 

0 . 4204 1 5 

42 

-0 . 02 96 07 

0.019078 

92 

0.006902 

0.372005 

4 3 

-0 . 026 425 

0 . 062426 

93 

0.007548 

0 . 3407 57 

44 

-0.024004 

0 .095400 

94 

0 .008408 

0 . 337226 

45 

-0 . 022 344 

0.117956 

95 

0.009775 

0.343233 

4 6 

-0 . 02 0640 

0 .1 40685 

96 

0.011507 

0.298582 

4 7 

-0.019004 

0.163511 

97 

0.013561 

0.242132 

48 

-0.017246 

0 . 186 925 

98 

0.015451 

0.209926 

4 9 

-0 . 0 1 5446 

0 .211437 

99 

0.017280 

0 . 179721 

50 

-0.013567 

0.237579 

100 

0.018999 

0 . 150880 


ORIGINAL PAGE IS 
OF POOR QUALITY 


SEG 

S 

BETA 

101 

0 . 020674 

0 . 1 22772 

102 

0 . 022342 

0.094780 

103 

0.023998 

0 . 06 6 99 1 

104 

0 .026420 

0 . 026 342 

105 

0 . 029605 

0.000000 

106 

0 . 032716 

0.000000 

107 

0 . 035768 

0.000000 

108 

0 . 038847 

0.000000 

109 

0 . 044215 

o . oooooo 

110 

0.051835 

0.000000 

i ; : 

0 . 059429 

0.000000 

112 

0 . 066997 

0.000000 

113 

0 .07 4547 

0.000000 

114 

0 . 082086 

0 .oooooo 

115 

0 . 089618 

o . oooooo 

116 

0 .097143 

0.000000 

117 

0.104666 

o oooooo 

118 

0.112187 

0.000000 

1 1 9 

0.119707 

0.000000 

120 

0 . 127227 

0 . 0 J 0 0 0 0 

121 
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HACA 0012 : EXAMPLE t : TIME* 120.000 SEC 


ICING CONDITION: 


* 


STATIC TEMPERATURE (C) 260.55 
STATIC PRESSURE (PA) 90744.00 
VELOCITY (M/3) 129.00 
LUC (G/n-*3) 0,50 
DROPLET DIAMETER (MICRONS) 20.00 


ICE ACCRETION DATA: 


STAGNATION POINT 74 

TRANSITION POINTS ( LOWER , UPPER ) 72 75 

ICING LIMITS ( LOWER, UPPER) 43 103 

NUMBER OF POINTS 147 

NUMBER OF SEGMENTS ADDED 4 


Figure 7.5: continued 
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Figure 7.5: continued 
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Figure 7.6: 
Example 1. 
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NACA 0012 : EXAMPLE 2 

ES26Y 
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Figure 7.7: Input data file for Example 2 
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Figure 7.10: y 0 vs. s data calculated in subroutine COLLEC for Example 2 
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Figure 7.1 It Icing parameter plots for the first timestep of Example 2. 
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Figure 7.11: Concluded 
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Figure 7.12: Locations of the multi; 
step of Example 2 




Figure 7.13: Locations of the multiple stagnation points in Example 2 after 
specifying new axes limits 
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Figure 7.15: yo vs. s data calculated in subrouting COLLEC for the second 
time step of Example 2 
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Figure 7.16: Icing parameter plots for the second timestep of Example 2. 
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Figure 7.16: Concluded 
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Figure 7.17: Comparison of experimental and calculated ice accretion 

shapes for Example 2. 
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Figure 7.18: Input data file for Example 3. 
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Figure 7.19: Icing parameter plots for the first timestep of Example 3. 
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Figure 7.19: Concluded. 
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Figure 7.20: Icing parameter plots for the second timestep of Example 3. 
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Figure 7.20: Concluded 
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Figure 7.21: Comparison of experimental and calculated ice accretion shape 
for example 3. 
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Figure 7.22: Droplet distribution specified for Example 4. 



Figure 7.23: Cumulative volume fraction vs. droplet diameter for the 
droplet distribution specified in Example 4. 
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ILIFT* 1 
IPARA* 1 
IFIRST* 3 
ISECND* 3 
IPVOR* 1 
IMCLT* 0 
CLT* 0.0 
ICHORD= 0 
OCl* 0.0 
IHD= 1 
ISOL* 0 
IPRINT* 0 
IFLLL* 1 
CEND 
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0.0318550 
0 .0199080 
0 .0000000 
CTRAJ1 
GEPS* 0 
DSHIFT* 
VEPS* 0 
0 
0 
1 
0 
1 
1 
1 


EXAMPLE 9 


0 . 9899999 
0 .8750000 
0.7250000 
0 . 5750000 
0 .9250000 
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0.0015000 
0 . 0050000 
0 . 0 3.0 0 0 0 0 
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0 .1750000 
0.3250000 
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-0.0179530 
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-0 . 097 9060 
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^ . 0600760 
■0 .0507320 
0.0356790 
0.0237260 
0.0086200 
0 .0051960 
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0 .0039620 
0.0069850 
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0 .0910380 
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0 . 0600990 
0.0596980 
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0.0000000 
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0.9800000 
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0.5500000 
0.9000000 
0 .2500000 
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0 .0950000 
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0 . 5000000 
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-0 . 0036110 
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-0.0997930 
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CEND 

CTRAJ2 
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CDIST 
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CICE 
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LWC* 0.750 
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PAMB* 90798.0 
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Figure 7.24: Input data file for Example 4. 
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Figure 7.25: Icing parameter plots for the first timestep of Example 4. 
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Figure 7.25: Continued. 
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Figure 7.25: Concluded. 
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Figure 7.26: Comparison of the calculated local collection efficiency vs. 
surface distance for a normal and monodispersed droplet distributions with 
mass median diameters of 20.0 microns. 
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Figure 7.27: Electrothermal heater on a NACA 0012 airfoil modeled in 
Example 5. 
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Figure 7.28: Input data file for Example 5. 
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Figure 7.29: Icing parameter plots for the first timestep of Example 5. 
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Figure 7.29: Concluded. 
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Chapter 8 

PROGRAMMING AIDS 

This section contains information regarding the coding of LEWICE to aid 
users who require a more advanced understanding of the program. Descrip- 
tions of the subroutines, common blocks, and work files are included, along 
with flow charts of various sections of the program. 

8.1 Descriptions of Subroutines 

The subroutines will not be discussed individually in this user’s manual. 
Instead, a description can be found in the program listing in COMMENT 
statements preceding most subroutines. These descriptions include the pur- 
pose of the subroutine, the input and output variables, and additional notes 
describing special features of the subroutine. 

8.2 Diagnostic and Error Messages 

In a program as complex as LEWICE, the programmer must provide mes- 
sages informing the user of any abnormalities in the calculations. The most 
common error messages are discussed in this section. They are grouped by 
the program module in which they will occur, i.e., potential flow calcula- 
tion, particle trajectory calculation, or thermodynamic and ice accretion 
calculation. 

When an error message or unusual situation occurs, the first step should 
be to check the description of the input parameters found in Section 5.0. 
The following descriptions of the error messages will assume that this step 
has already been completed, and will concentrate on how the error may 
relate to certain aspects of the icing condition or geometry being evaluated. 

8.3 Potential Flow Calculations 

All of the lines of data input to the potential flow code are identified by 
a card type number found in the most right-hand column in each line of 



data. When the data line is read, this value is checked against the value 
of the card type that was anticipated. If these numbers do not match, an 
error message is printed and the program terminated. This error condition 
is usually caused by having input values in the wrong columns on a line of 
data. 

Another error occurs occasionally when stagnation points are calculated 
behind (downstream) of the large horns of a glaze ice accretion. The error 
will usually be characterized by “division by zero" messages. The situation 
cannot be corrected by forming a pseudo-surface as described in Section 7.2 
and Appendix C. However, the case should be re-run with a larger timestep 
which will make the predicted glaze ice accretion somewhat smoother and 
may allow the calculations to continues. 

Another error condition can occur when multiple stagnation points are 
calculated which cannot be removed by the formation of a pseudo-surface, 
and a stagnation point has been manually selected. In this case, there will 
be at least two locations where the local velocity is 0.0 m/s. When the 
compressible, dimensional surface velocity is calculated from the incom- 
pressible, non-dimensional values (subroutine VEDGE), division by zero 
errors can again occur. It is best to force the calculations to continue 
through this error, if possible, because the locations of the error may be 
downstream of the impingement region and, therefore, will not affect the 
ice accretion shape. If the calculations cannot be continued, re-run the case 
with a larger timestep. 

8.4 Particle Trajectory Calculations 

Most of the errors that occur in the particle trajectory calculations result 
from inaccuracies in the flow field or errors in input data. Errors caused 
by improper input parameters will generally be accompained by diagnostic 
messages, which will be discussed below. 

Subroutine RANGE will allow a total of 30 trajectories to be calculated 
while searching for a trajectory that passes above and below the body. 
If two such trajectories are not identified, the following message will be 
printed and the program stopped: 

30 Trajectories are calculated in RANGE. Run aborted. 
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The most probable cause of the error is that the values of YOMAX and 
YOMIN have been input incorrectly. If the values are correct, review the 
body coordinates that are being used by the program. 

A run can be terminated in a similar manner in subroutine IMPLIM. In 
this case, if the number of trajectories required to identify the upper and 
lower impingement limits exceed NSEAR, the run is terminated. The value 
of NSEAR is input by the user through namelist TRAJ1 and is normally 
set equal to 50. If the limiting value of NSEAR is reached, the value of 
YOLIM may be unnecessarily small or there may be a problem with the 
way the program has read or interpreted the body coordinates. 

In subroutine COLLEC, the program will be terminated if the value of 
NPL is greater than 100. The value of NPL is input by the user through 
namelist TRAJ1, and typical values are between 10 and 20. If this error 
occurs, verify that the value of NPL is in the proper column in the input 
file. 

Occasionally, errors occur in the flow field near the surface of the body, 
especially for convoluted glaze ice accretions. In these cases, very large 
local velocities are calculated which cause exponent overflows or negative 
arguments in subroutine ABFORM. If these errors occur, try to force the 
calculations to continue through this error because, if a bad impingement 
point is calculated as a result of the error, it can removed before calcu- 
lating the local collection efficiency. Unfortunately, once this error has 
been encountered in one timestep, it will probably occur in all subsequent 
timesteps. 

8.5 Thermodynamic and Ice Accretion Calculations 

If the program completes the flow field and particle trajectory calculations 
with no errors, there generally will be no additional errors in the ther- 
modynamic calculations. There are temperature limits in the subroutines 
used to calculate the pressure of water vapor over liquid water and over 
ice. These calculations are performed in subroutines PVW and PVI, re- 
spectively. The temperature ranges, shown below, are sufficient for any 
anticipated application of the code. 

Vapor Pressure over Liquid Water: 223.15 K < T < 323.15 K 
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Vapor Pressure over Ice : 


213.15 K < T < 273.15 K 


8.6 COMMON Blocks and Work Files 

In a program of this size, much of the information must be transferee! 
between subroutines through COMMON statements. Table 8.1 lists each 
COMMON block in the program, the general purpose of the variables in 
the COMMON block, and the subroutines in which each block is found. 
No open COMMON statements are used in LEWICE. 

In addition to COMMON blocks, much information is passed between 
subroutines through temporary work files. The work files and the purpose 
of each are shown in Table 8.2. Many of these files are used in the potential 
flow calculations (subroutine S24Y), and will not be described in detail 
because little work was done to modify the potential flow calculations. 

8.7 Size of the Code 

Because LEWICE combines three complex computer codes into a single 
computational algorithm, the code requires a substantial amount of com- 
puter memory for both the source code and operation. On the IBM 370 
computer used at NASA Lewis, the source code itself requires approxi- 
mately 1400 K of memory. The input file requires only 4 K, but 500 K 
should be allowed for the output file, especially if a droplet distribution has 
been specified. Also, a restart file is generated by the program which will 
be the same size as the input file. The work files used in the code require 
approximately 150 K of memory. They are erased at the termination of the 
run and, therefore, could be placed in a temporary work area. 
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Table 8.1 Common Blocks Used in LEWICE 


COMMON Block 

Description 

LABELS 

Contents: 

Variables: 

Subroutines: 

Parameters to be placed on plots 

IDR, VINF, LWCKG, TAMB, PAMB, RH, 

MAIN, BORDER 

CNTL 

Contents: 

Variable: 

Subroutine: 

Control parameters used in the plotting 
and particle trajectory routines 
LOPT, LEQM, IP LOT, LCMB, LCMP 
MAIN, TRAJ, IMPLIM, COLLEC, INTIG, 
RANGE, MODE, READIN, VELCTY, COMB2I) 

PSEUDO 

Contents: 

Variables: 

Subroutines: 

x-coordinates of the pseudo surface that have 
been generated 

IPSURF XPS 

MAIN, STAG, VEDGE, PSURF, NEW45, 
MODE, READIN 

ICECOM 

Contents: 

Variables: 

Variables used in the thermodynamic and 

ice accretion subroutines 

X, Y, SEGLEN, SEGLIN, VE, TE, PE, 


RA, XK, BETA, TSURF, FFRAC, HTC, RI 
QCOND, MDOTC, MDOTE, MDOTRI, MDOTTI, 
MDOTT, DICE, VINF, TAMB, PAMB, RH, 
DPMM, LWCKG, CPA, CPI, LV, LF, VISC, 

PI, NPTS, NTHI, NTLOW, ISTAG 
Subroutines: CNSTS, STAG, PSURF, NWPTS, SEGSEC, 

ICE, CDYLYR, EBAL, COMPF, COMPT, 
NWFOIL, OUTPUT, PLOTD 


150 


Table 8.1 Continued 


COMMON Block Description 

DROP Contents: Droplet distribution data 

Variables: NSI, FLWC, IRS, DPD, EPT 

Subroutines: PLOTD, TRAJ, EFFICY 

BFLAG Contents: Potential flow data 

Variables: IDB, INL, IFL, NL, LIFT, IBMF, ISAV1, 

ISAV2, ISAV3, BTITLE, IBT, IBST, 

IBTOT, NELTOT, ITRB, INMB, CHORDB, 

IBD, LIFTOT, IPRB, IFST, ISEC, 

FTITLE, IPVR 

Subroutines: S24Y, MAINl, MAIN3, ASSEMB, ELFORM, 

FLOWS, MAFORM, PRNTEL, VXYOFF 

COMBO Contents: Potential flow data 

Variables: CCL, INCLT, CLT, ALPHA, SUMDS, TLU, 

IND, ALPHAO, CNU, SMDSWF, MIO 
Subroutines: S24Y, MAINl, MAIN3, ASSEMB, COMBO, 

FLOWS, MAFORM, OFFBOD, VXYOFF 

FILEID Contents: Potential flow data 

Variables: IFILE1, IFILE2, IFILE3, IFILE4, 

IFILE5, IFILE6, IFILE7, IFILE8, 

IFILE9, IFIL10, IFIL11, IFIL12, 

IFIL13, IFIL14, IFIL15, IFIL16, 

IFIL17, IFIL18, IFIL19, IFIL20 
Subroutines: S24Y, REWYND, FILRS, MAINl, MAIN3, 

SOLVE, ASSEMB, COMBO, ELFORM, FLOWS, 
MAFORM, VPROFF, VXYOFF 

MDATA Contents: Potential flow data 

Variables: ISOL, IOFF, NONU, MBNU, IPRINT, MORE, M 

Subroutines: S24Y, MAIN1M MAIN3, ASSEMB 
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Table 8.1 Continued 


COMMON Block 

Description 

ROTAT 

Contents: 

Variables: 

Subroutines: 

Potential flow data 
NROT, ROTRAD 

S24Y, ASSEMB, FLOWS, MAFORM 

ELDATA 

Contents: 

Variables: 

Subroutines: 

Potential flow data 

XO, YO, DS, SA, CA, CURV, DL 

MAIN1, ASSEMB, ELFORM, MAFORM 

GCOEFS 

Contents: 

Variables: 

Subroutines: 

Potential flow data 

CD, CF, CG, Cl, WF 

MAINl, ASSEMB, ELFORM, MAFORM 

COM 

Contents: 

Variables: 

Subroutines: 

Potential flow data 
IFILL 

MAINl, MAIN3, FLOWS, VPROFF, VXYOFF 

SPACER 

Contents: 

Variables: 

Subroutines: 

Potential flow data 

WKAREA 

SOLVE 

SIGMAS 

Contents: 

Variables: 

Subroutines: 

Potential flow data 
CSIG, CK 

COMBO, FLOWS, VXOFF 

GEOMD 

Contents: 

Variables: 

Subroutines: 

Potential flow data 
X, Y, XSAVE, YSAVE 
ELFORM, FLOWS, BTITLE 

GCFF 

Contents: 

Variables: 

Subroutines: 

Potential flow data 
CD, CF, CG, Cl, WF 
FLOWS, VXYOFF 
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Table 8.1 Continued 


COMMON Block 

Description 

ELDD 

Contents: 

Potential flow data 


Variables: 

X, Y, DS, SA, CA, CURV, DL 


Subroutines: 

FLOWS, VXYOFF, READIN, VELCTY 

OVER 

Contents: 

Potential flow data 


Variables: 

VI Y, V2Y, V3Y, V4Y, V5Y, V6Y, V3X, 
V4X, V5X, V6X 


Subroutines: 

VPROFF 

BODY 

Contents: 

Coordinates describing the body 
geometry 


Variables: 

NPTS, XSH, YSH, X, Y 


Subroutines: 

TRAJ, MODE, READIN, PLTRAJ 

BOUND 

Contents: 

Coordinates of the most forward and 
aft points of the body geometry and 
boundaries for the particle trajectory 
calculations 


Variables: 

XFRNT, YFRNT, XREAR, YREAR, XSTOP, 
YLO, YUP 


Subroutines: 

TRAJ, RELEAS, IMPLIM, INTIG, RANGE, 
MODE, READIN 

IMP 

Contents: 

Particle impingment data for use in 
the collection efficiency calculation 


Variables: 

S, SW, YO, IS 


Subroutines: 

TRAJ, IMPLIM, COLLEC, ORDER, INTIG, 
RANGE, MODE, READIN 

INIT 

Contents: 

Initial conditions for particle 
release 


Variables: 

XIN, YIN, DPM, RL, PITDOT, PIT, 
VINF, VXP, VYP, AOAR 
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Table 8.1 Continued 


COMMON Block 

Description 


Subroutines: 

TRAJ, RELEAS, IMPLIM, COLLEC, INTIG, 
RANGE, DIFFUN, MODE, COMB2D 

DIFF 

Contents: 

Variables used in the particle 
trajectory equations 


Variables: 

Q, AMASS, G, YYI, VISC, CF 


Subroutines: 

TRAJ, DIFFUN 

STEP 

Contents: 

Time step and error criteria used 
in the integration of the particle 
trajectory equations 


Variables: 

TIMSTP, EPS 


Subroutines: 

TRAJ, IMPLIM, COLLEC, RANGE 

TP15 

Contents: 

Variables used to calculate the 
combination solution of velocity 


Variables: 

ABCD, ATOTAL, VC15, RSORTC, NX15 


Subroutines: 

READIN, COMB2D 

NUM 

Contents: 

Variables used to calculate the 
flow field velocities 


Variables: 

N, M, IND, IBTOT, LIFTOT 


Subroutines: 

READIN, VELCTY 

GCF 

Contents: 

Parameters used to calculate the 
flow field velocities 


Variables: 

CD, CF, CG, Cl, WF, CSIG, CK, SIG 


Subroutines: 

READIN, VELCTY 

FLG 

Contents: 

Parameters used to calculate the flow 
field velocities 


Subroutines: 

READIN, VELCTY 
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Table 8.1 Continued 


COMMON Block 
VEX Contents: 

Subroutines: 


Description 

x- and y-components of the flow 
field velocity 
VELCTY, COMB2D 


Table 8.2 Output and Work Files Used in LEWICE 


File 

Number Type Description 


1,2,4, 

7-18, 

21 


These are work files used in the potential 
flow calculations. Since no modifications 
were made to S24Y, the contents of these files 
have not been examined in detail. 


3 Binary 


5,6 Character 

20 Character 


22 Character 


Unit 3 contains the yo vs. s points for each 
droplet diameter. This file is read in 
subroutines EFFICY and PLOTD. 

These are the read/write defaults to display 
information on the screen. 

This file contains the airfoil coordinates for 
all timesteps. The file is read in subroutine 
PLOTD to plot the airfoil. 

Unit 22 contains the airfoil segment distances 
used in the collection efficiency 
calculations. 
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Table 8.2 Continued 


File 

Number 

24 

30 

35 

36 


45 

56 


Type 

Character 


Binary 


Character 


Character 


Character 


Character 


Description 

This is an internal input file for subroutine 
TRAJ. It is created from the input on 
unit 35. There is no updating of unit 24 for 
the second or subsequent timesteps. 

This is the main working file for the 
program. It contains the current values of 
the body coordinates, surface distances, 
edge velocities, pressures, etc. for each 
segment. 

This is the primary input file prepared by 
the user. It is used to set up all of the 
other internal input files. 

This is the restart input file produced by 
the program after each timestep. It is 
identical to the input file on unit 35, 
except that the body coordinates are of the 
current ice shape. The icing time 
corresponding to the ice shape coordinates 
is also included in namelist ICE. 

This is the internal input file for the 
potential flow code (S24Y). It is formed from 
the data input on unit 35. 

This is the primary output file for LEWICE. 
It contains the printed output from all 
portions of the code. 
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8.8 Execution Times 


When an ice shape is predicted, the largest portion of the execution time is 
spent on the particle trajectory calculations. The CPU time required to run 
the example cases found in Section 7.0 on the IBM 370 computer are given 
in Table 8.3. Also indicated are the total number of particle trajectories 
calculated in each, and the average CPU time per trajectory. 

Table 8.3 CPU Time for Examples 1-5 


Example 

CPU Time 

Number of 
Trajectories 

CPU-sec 

Number 

(sec) 

Calculated 

Trajectory 

1 

362.52 

67 

5.41 

2 

350.80 

71 

4.90 

3 

395.01 

68 

5.81 

4 

719.30 

205 

3.51 

5 

185.50 

44 

4.22 


With these results, the computational time required for a specific icing 
condition can be estimated if the speed of another type of computer is 
known compared to an IBM 370. 

8.9 Getting LEWICE Operational 

When LEWICE is transferred to a system other than that at NASA Lewis, 
it must first be compiled. A FORTRAN 77 compiler has been used to 
compile the version used at NASA Lewis. The graphics commands used 
in LEWICE are from a system that is unique to NASA Lewis and will 
not be directly applicable to any other system. The name of this graphics 
package is GRAPH3D, but most of the commands in LEWICE are from 
GRAPH2D, the predecessor of GRAPH3D. Calls to the graphics package 
are located in subroutines STAG, PSURF, PLOTD, BORDER, INTIG, 
PLTRAJ, and EFFICY. The graphics in LEWICE use only the simplest 
commands, and, therefore, modification of the commands should not be 
difficult; unfortunately, it will probably be time-consuming. 
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Appendix A 

DERIVATION OF THE ICING ENERGY EQUATION 

In this appendix, the thermodynamic processes that take place during 
the formation of an ice accretion are identified and expressed mathemati- 
cally to form the icing energy equation. This derivation is similar to that 
presented in Reference A-l, except that the equation is applicable to an 
arbitrary control volume instead of a control volume situated on the stag- 
nation line. 

A.l Definition of the Control Volume 

The control vol um e to be analyzed is located on the surface of the body 
and extends from outside the boundary layer to the surface of the body, 
as shown in Figure A-l. It encloses a distance along the external surface, 
and, for dimensional completeness, extends one unit length in the spanwise 
direction (into the page). The lower boundary of the control volume is 
initially on the surface of the clean geometry, and moves outward with the 
surface as the ice accretes. Therefore, the control volume is always situated 
on either the clean or iced surface, and any accumulated ice is considered 
to leave the control volume through the lower boundary. This definition is 
important in conduction heat transfer, discussed later in this appendix. 

A. 2 Mass Balance on an Icing Surface 

An evaluation of all mass entering and leaving the control volume is 
shown in Figure A-2. A mass balance equation can be formed from these 
terms as shown below. 


m c + m r . n -m e - m rnl = m, (A — I) 

At the stagnation point, there will be no water inflow along the surface 
so therefore, m r . n = 0.0. 

Since the freezing fraction is defined as the proportion of the total mass 
of liquid entering the control volume that freezes in the control volume, it 
can be expressed by the following equation: 


/ = 


rh, 

rh e + m r . n 


(A-2) 
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Substituting Equation (A-l) into Equation (A-2), the water flow out of 
the control volume can be expressed as 


= (1 - f)(rh e + fhrj - rh, (A- 3) 

A. 3 Energy Balance on an Icing Surface 

The same control volume concept is used to formulate the energy bal- 
ance on the icing surface. The First Law of Thermodynamics for a control 
volume can be expressed as: energy inflow rate = energy outflow rate + 
energy storage rate. 

The modes of energy transfer, illustrated in Figure A-3, are as follows: 


Mode of Energy Transfer 

1. Impinging Water 

2. Water Flow Into Control Volume 

3. Evaporation 

4. Water Flow Out of Control Volume 

5. Ice Accumulation Within Control Volume 

6. Convection 

7. Conduction through the Skin 


Energy Flow Rate 

m c i wJ 

rhii i<tur 

qeAs 

qkAs 


Using the convention that energy flow into the control volume is positive, 
the terms can be summed to yield the general form of the energy equation: 

rh e tvi.T “I - di r<ft t‘w,«ur(t— i) = t V|(ur *4* m rom , t w tur 

+ rhi i iitur + q e As + q k As ( A - 4) 

The evaluation of the terms of the energy equation has been made by 
various authors, most notably by Sogin (Reference A-2), Lowzowski (Refer- 
ence A-3), and Cansdale and Gent (Reference A-4). The following sections 
will evaluate each of the terms of Equation (A-4). 
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A. 4 Impinging Water 

Since droplets are essentially brought to rest when they strike an object, 
it is appropriate to use the stagnation enthalpy defined as 

T = (T. - 273.15) + 'f (A- 5) 

The arbitrary reference for zero enthalpy used in this study is water at 
273.15 K. Substituting Equation (A-5), the energy flow rate of the imping- 
ing water becomes 

m e i.,r = m’ [c,„, (T. - 273.15) + ^ ] As - 6) 

where 

m” = LWC{V , „)/? {A - 7) 

A. 5 Water Flow Into the Control Volume 

The water flowing into the control volume will be at the surface tem- 
perature of the preceding control volume. The enthalpy can therefore be 
expressed as 

*w p »ur(»-l) = C P«, 1 ) (2»«r(»-I) “ 273.15) (A — 8) 

where (i-1) denotes that the specific heat and temperature are evaluated at 
the conditions of the proceding control volume. The ronback water energy 
flow rate into the control volume can therefore be expressed as 

rhr in V, ur(i-i) = "»r„, Cp., (T^i-i) - 273.15^ (A - 9) 
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A. 6 Evaporation 

The rate of energy transfer from the surface because of evaporation is 
given by 


r — (P*w 273.15) + ■i'v j As (A 10) 

where m" is the evaporative mass transfer flux and L v is the latent heat of 
vaporization. 

The mass transfer rate is analogous to the convective heat transfer rate 
and can be written as 


rh" = g AB 


(A - 11) 


where g is the mass transfer coefficient and A B is the evaporative driving 
potential. The mass transfer coefficient, g, can be evaluated using the 
analogy to heat transfer given by the following equation found in Reference 
A-5: 



(Pr \ M1 

\Sc) 


[A - 12) 


In Equation (A-12), Pr is the Prandtl number and Sc is the Schmidt num- 
ber. 

The mass transfer driving potential is analogous to the temperature 
difference in the convective heat transfer equation. In the case of evapora- 
tion, the driving potential is a vapor concentration difference instead of a 
temperature difference. The equation used in this study is similar to that 
derived by Sogin (Reference A-2), and given as 


Pv,mr /Pfur r h P\i,$/ P$ 

A B = T ^ (A - 13) 

(l/. 622 jP e /T e - Pv.tur/Ttur 
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This term accounts for compressibility effects, as does the term derived 
by Cansdale in Reference A-4. 

When the water droplets freeze on impact (f = 1.0), there is no liquid 
water on the surface to be evaporated; however, water vapor can still leave 
the surface through sublimation. In this case, Equations (A-10), (A-ll), 
(A-12), and (A-13) are still used to determine the rate of energy transfer 
from the surface, except that the latent heat of sublimation, L,, is used in 
Equation (A-10) instead of the latent heat of vaporization, L v . 

A. 7 Water Flow Out of the Control Volume 

The water flowing out of the control volume will be at the surface tem- 
perature of the control volume, allowing the enthalpy to be expressed as 

*m,»ur(») = c p„, — (T,ur(.) _ 273.15) (A — 14) 

where (i) denotes that the specific heat is to be evaluated at the surface 
temperature of the control volume being analyzed. Using Equation (A-3), 
the runback water energy flow rate can be expressed as 




.) = (!-/)(" 


me +m ri J— m. 


U P* t .*r(0 


(r w( ,)-273.15) (A -15) 


A.8 Ice Accumulation Leaving the Control Volume 

As previously discussed, the control volume remains on the surface of the 
geometry as the ice accumulates within the control volume. 

From the definition of the freezing fraction, Equation (A-2), the freezing 
rate is 


rhi = f (m e + m r .J 

The enthalpy of ice referenced to water at 273.15 K is 

= c p;,.»r {T$ur ~ 273.15) — Lf 


(A - 16) 
(A - 17) 


164 



Combining Equations (A-16) and (A-17), the rate of energy leaving the 
control volume in the accumulated ice can be expressed as 


*hi — f {pic + *hf;„ ) 273.15) — Zf/J (A — 18) 

A.9 Net Convective Heat Flux 

The local value of the aerodynamically induced heat flow to or from 
the outer boundary of the control volume is determined by the convective 
cooling and kinetic heating of the surface. The net convective heat flow 
can therefore be defined by the following equation: 

9c As = h e (T tur T a% ,) As (A — 19) 

In this equation, T tur is the surface temperature and T aw is the adiabatic 
wall temperature. The adiabatic wall temperature is given by 

V 2 

T aw = T e + r e — (A - 20) 
2 

where T e and V t are the temperature and velocity at the edge of the bound- 
ary layer, respectively, and r c is the recovery factor. The local temperature 
is calculated from the pressure coefficients calculated by the potential flow 
code using the isentropic relationships. Substituting Equation (A-20) into 
Equation (A- 19), the expression for the convective heat flow rate becomes 

9c As= h e {T. w -T e - r ^-)As (A -21) 

i c Pm 

The heat transfer coefficient is calculated using the integral boundary 
method described in Appendix B. 

A. 10 Conduction From the Airfoil Surface 

When the cloud is first encountered, a temperature difference will exist 
between the wetted surface and the inner structure of the airfoil. Prior to 
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entering the cloud, this inner structure is assumed to be at an equilibrium 
temperature. The evaluation of the resulting conductive heat flow rate is 
dependent on knowing the thermal conductivity and detailed geometry of 
the inner structure of the airfoil. 

After a layer of ice forms on an unheated surface, the temperature of 
the skin should again reach an equilibrium temperature. Since ice is an 
insulator, any heat transfer through the skin will not affect the growth of 
the ice accretion at the air/ice interface. 

In a thermal deicing system, an ice layer is allowed to form and heat is 
then applied at the ice/skin interface. The effect of this heat is to melt the 
ice attached directly to the surface, thereby allowing the ice to shed due 
to the aerodynamic forces acting on the ice accretion. Currently, LEWICE 
is not capable of analyzing this deicing phenomenon because a thermal 
analysis would be required not only at the air/ice interface but also at the 
ice/surface interface. A complete description of the current capability to 
model a thermal deicing system is given in Reference A-6. 

A thermal anti-icing system differs from a deicing system in that suffi- 
cient heat is applied to prevent any ice from forming. The formulation of 
the icing energy equation in LEWICE should be applied to thermal anti- 
icing systems only to determine the minimum heating requirements to keep 
ice from forming on the surface. The heat flux from the skin is specified as 
a function of s and is assumed to be from the inner control volume bound- 
ary. On the first timestep, this boundary is the uniced surface. If the heat 
input is not sufficient to keep ice from forming, the heat input on the second 
timestep would be assumed to come from the iced surface and incorrectly 
neglects the thermal conductivity of ice and the effect of the varying ice 
thickness on the surface. An example of the application of the code to a 
thermal anti-icing system is discussed in Section 7.5. 
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The energy terms can now be summed to form the complete energy 
balance for an icing surface used in LEWICE. This equation is as follows: 


™c [cp.„ {T, - 273.15) + ^] + 

(7’«ir(i-i) — 273.15) j + qt As = 
m e [c Pmnr (Tsur - 273.15) + L v j + 

[(1 - f){rh e + m r J - m.J c Pwnr {T^ - 273.15)+ 
f(rh e — m rt(> ) {T,ur — 273.15) — Lf J + 

he [t. w - T t - As (A - 22) 

The solution of Equation (A-22) is presented in Section 4.3.2. 
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Figure A-3: Energy balance for a control volume. 


APPENDIX B 


INTEGRAL BOUNDARY LAYER METHOD ON AN ICED 
SURFACE 

As described in the text, the evaluation of the boundary layer character- 
istics using a complete Navier-Stokes solver would be too time-consuming 
to be practical in an interactive program such as LEWICE. An integral 
boundary layer solution was therefore developed to account for surface 
roughness and variable velocity in the calculation of the convective heat 
transfer coefficient. The original algorithm was developed by UDRI, but 
was modified in the current version of LEWICE. 

B.l Characterization of the Ice Surface Roughness 

The integral boundary layer method to be applied in LEWICE is re- 
quired not only to determine the convective heat transfer coefficient on the 
iced surface in the laminar and turbulent regions, but also to determine the 
location of the transition from laminar to turbulent flow (transition point). 
The transition to turbulent flow is caused primarily by the surface rough- 
ness of the iced surface which acts to trip the boundary layer from laminar 
to turbulent flow. Also, once turbulent, the surface roughness elements 
enhance the convective heat transfer coefficient by 1) increasing the skin 
friction coefficient and 2) increasing the effective surface area from which 
heat transfer takes place. 

A representative surface of an actual ice accretion is shown in profile in 
Figure B-l. The size and shape of the surface roughness elements found on 
typical ice accretions are functions of the atmospheric and meteorological 
conditions at which the ice is formed. Unfortunately, there is insufficient 
data to characterize the size and shape of the surface roughness as a func- 
tion of those conditions. Also, the complex roughness patterns found on 
typical glaze ice accretions are beyond the analysis capability of an integral 
boundary layer method. Therefore, an alternate method to characterize 
the surface roughness is required. 

Analysis of rough surfaces has been performed by numerous investiga- 
tors, one of the earliest and most well-known being Nikuradse (Reference 
B-l). These experiments dealt primarily with turbulent flow in pipes that 
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were artificially roughened with uniform grains of sand. Schlichting (Ref- 
erence B-2) introduced the concept of equivalent sand-grain roughness as 
a means of characterizing other types of roughness elements by referring 
to the equivalent net effect produced by Nikuradse’s experiments. In this 
application, the irregular roughness elements of an iced surface are repre- 
sented by a value of equivalent sand-grain roughness height, k,, as shown 
in Figure B-2. This value is specified by the user and is constant during the 
calculation of a given icing condition. The value of k, should be changed as 
a function of the atmospheric and meteorological parameters of the icing 
condition using the guidelines discussed in the text. 

B.2 Calculation of the Boundary Layer Characteristics 

Determination of the Transitio n Location 

The evaluation of the boundary layer is begun at the stagnation point 
(s = 0.0). The calculations are initiated by determining the location of the 
transition for laminar to turbulent flow, i.e., the transition point. The crite- 
ria for transition, developed by Von Doenhoff (Reference B-3), assumes that 
the flow will become turbulent when the local roughness Reynolds number 
is greater than 600. The empirical criteria, based on data obtained using 
sand-grain type roughness elements, are therefore given by the following 
equation 


R = Yl*± >600 {B - 1) 

* V 

where k, is the equivalent sand-grain roughness height representing the 
actual ice surface roughness. The velocity at y = k , , designated V*, (Figure 
B-2) must be determined to evaluate the local roughness Reynolds number. 
The following is the method used to determine V k , found in Equation (B-l). 

Assume that the laminar velocity profile can be represented by a 4th 
order polynomial of the form 

= «. + ..( f) + «,(«)* + «,(§)• + «4 r (b - 2 ) 
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where 6 is the boundary layer thickness, and V t and y are defined in Figure 
B-2. This assumption is known as the Pohlhausen approximation (Refer- 
ence B-4, pp. 310-311). By applying the following boundary conditions at 

y = 0: 


d 2 V 1 dp 

V dy 2 p da 

dny n 
lim -r— = 0 
y-K» oy n 

it can be shown that the expression for V /V e resulting from Equation (B-2) 
is 


i 

ii 

(£— 3a) 

n = 1, 2, 3, ... 

(B - 36) 


^ = [2(y/<) - Av/6) Z + 2(J </«)‘] + l/« A (?/<) " y/ ( ) a ( fl - 4 ) 

V t 

where 

A = «V„^! (B - 5) 

da 

If we let y = k t in Equation (B-4), the velocity at y = k t can be written 
as follows: 


y = [*(*,/«) - 2 (*./«)* + (A./*) 4 ] + 1/6 A (*./f)(l - *./«)» (fl - 6) 

To apply this equation, the boundary layer thickness, 6, must be evaluated. 

The laminar momentum thickness, 0j, is defined in Reference B-4, p. 
244, as 



Substituting Equation (B-4) into Equation (B-7) and integrating from 
y = 0 to y = 6, it cam be shown that 6 is related to the laminar momentum 
thickness by the following approximation: 

6 « 8.5 0i ( B - 8) 

From the integral momentum equation, the laminar momentum thick- 
ness can be evaluated using Thwaites formula (Reference B-4, p. 315) which 
is given by the following equation: 

* ?= i fL‘ v ’ d ‘ {B ~ 9) 

If it is assumed that the velocity at the edge of the boundary layer, V e , 
is the surface velocity calculated by the potential flow equations, Equation 
(B-9) can be numerically evaluated to determine for each segment on 
the body. Equations (B-6) and (B-8) can then be applied to determine V* 
as a function of s, and, therefore, Equation (B-l) can be evaluated along 
the surface to determine whether the flow has transitioned from laminar to 
turbulent flow. 

Laminar Convective Heat Tr ansfer Coefficient 

If the boundary layer is found to be laminar at a surface distance, s, 
the convective heat transfer coefficient for flow over a constant temperature 
body of arbitrary shape is calculated using an equation developed by Smith 
and Spaulding and described in Reference B-4, pp. 327-329. This equation 
is given as 


hi{s) = .296 





where A is the thermal conductivity of air. 
Turbulent Convective Heat Transfe r Coefficient 


(B ~ 10 ) 
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If, upon evaluating Equation (B-l), it is found that the boundary layer 
is turbulent, Equations (B-7)-(B-10) are not applicable, and an alternate 
method to determine the convective heat transfer coefficient must be de- 
veloped. Using a technique outlined in Reference B-5, an overall Stanton 
number can be developed using the thermal and momentum laws of the 
wall for fully rough flow. This equation is given in Reference B-5, p. 132 
as 


St = 


Sill 

Pr t + \Jcf/2 (1 /Stk) 


(B - 11 ) 


The terms of this equation that must be evaluated are the skin friction 
coefficient, c/, and the roughness Stanton number, Sf*. The experimental 
data for air (Reference B-6) suggest that the turbulent Prandtl number, 
Pr t , is approximately constant and equal to 0.9. 

Assuming for now that the values of cj and S'/* are known, the turbu- 
lent heat transfer coefficient is calculated using Equation (B-ll) and the 
definition of the Stanton number. Therefore, 


h t (s) = StpV e c p = 


cj[ 2 

,Pr t + Stu). 


pV t c p 


(B - 12 ) 


The expressions for the skin friction coefficient and the roughness Stan- 
ton number are now developed. 

Skin Friction Coefficient 


If the boundary layer has been found to be turbulent, i.e., the roughness 
Reynolds number is greater than 600.0, the surface can be considered to 
be fully rough (Reference B-5, p. 186). A basic characteristic of the fully 
rough surface is that the skin friction coefficient, c/, is independent of the 
Reynolds number. In this region, the pressure drag on individual roughness 
elements dominates and viscosity is no longer a significant variable. With 
this assumption, an expression for the skin friction coefficient can be devel- 
oped from the momentum law of the wall for fully rough flow (Reference 
B-5, pp. 186-188) as follows 
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(B - 13) 


£/ 

2 


.41 


1 * 


^ /n (86i^ + 2.568) 


where 0 t is the turbulent momentum thickness and k, is the equivalent sand- 
grain roughness height. (This equation was derived in a manner similar 
to Equation 10-45 in Reference B-5 except that Equation 10-43 was not 
simplified.) 

The turbulent momentum thickness is evaluated using the momentum 
integral equation in a manner similar to that for the laminar momentum 
thickness. The equation for 0 t is given in Reference B-5, p. 175 as 


0 t (s) 


'.0156 

v* M 



(B - 14) 


Since the turbulent boundary layer is proceeded by a laminar boundary 
layer, the numerical integration of Equation (B- 14) is begun at s = s tr 
instead of s = 0.0. The laminar momentum thickness already existing at 
s = s tr must then be added. Equation (B-14) can therefore be written as 




.0156 r* 
V* 11 J Sir 


1.8 


V c 3M ds 


+ 0l(str ) 


{B - 15) 


where 0[ is evaluated using Equation (B-9) and V t is the surface velocity 
calculated by the potential flow code. 

Roughness Stanton Number 

The roughness Stanton number is developed from the thermal law of 
the wall for fully rough flow, and is given in Reference B-5, p. 231 by an 
equation of the form 




(B - 16) 
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where V r is the shear velocity. In Equation (B-16), C is a constant that 
must be determined from experimental data and is a function of the type 
of roughness. Data for a rough surface composed of closely-packed spheres 
yields a value of C = 1.0 (Reference B-6). Setting the Prandtl number, 
Pr, equal to 0.72 and substituting the values for C and Pr into Equation 
(B-16), the expression for St* becomes 

St* = 1.16p£) (£-17) 

The shear velocity, u r , is evaluated using the equation from Reference B-5, 
p. 187. 

£ = (-B - 18) 

In this equation, cj is determined from Equation (B-13). 

B.3 Method of Solution 

To correctly apply the equations previously discussed, the integration 
procedure should be identified. As discussed in the text, the geometry is 
represented by a set of Cartesian coordinates (nodes) connected by straight 
line segments, as shown in Figure B-3. The stagnation point is designated 
s = 0.0, and will always fall on a nodel 1 point. The s values used in the 
integration correspond to the midpoint of each of the body segments. The 
calculation of the boundary layer characteristics is begun by evaluating 
Equations (B-6), (B-8), and (B-9) for s = 0.0 to s = s x . The transition 
criteria, Equation (B-l), is then applied, and, if the flow is found to be 
turbulent at this point, it is assumed to be turbulent at each segment 
downstream. Equations (B-13), (B-15), (B-17), (B-18), (B-ll), and (B-12) 
are applied to calculate h c if the boundary layer is turbulent. If the flow is 
laminar, Equation B-10 is used to calculate h e and the turbulence criteria 
is checked at the next segment. This procedure is then repeated for the 
lower surface. 
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B.4 Comparison with Experimental Data 

The results of the integral boundary layer method described in the pre- 
vious sections have been compared to experimental convective heat transfer 
data collected by Achenbach (Reference B-7). The current method was also 
compared to an integral boundary layer method developed by Makkonen 
(Reference B-8). 

An excellent discussion of the analytical method developed by Makko- 
nen and the comparisons to experimental data can be found in Reference 
B-8. The work by Makkonen also includes a discussion about the applica- 
bility of an integral boundary layer method to the calculation of convective 
heat transfer characteristics over an ice accretion shape. Therefore, a de- 
scription of the experimental data and detailed analysis of the results will 
not be presented in this appendix. Instead, the results of the integral 
boundary layer method described in this appendix will only be compared 
to the results in Reference B-8 and the general trends identified. 

The experimental measurements were made on a 15-cm diameter cylin- 
der roughened with grains of sand at various Reynolds numbers. The rough- 
ness element height, k, is 0.9 mm and the equivalent sand grain roughness 
height, k t , is 1.35 mm. The method developed by Makkonen uses both 
the maximum probable roughness height and the equivalent sand grain 
roughness height. The method described in this appendix uses only the 
equivalent sand grain roughness. Figure B-4 shows a comparison of the 
analytical results of LEWICE (with k t = 1.35 mm) and Makkonen to the 
experimental results. At a Reynolds number of 4.8 x 10 4 , the method of 
LEWICE shows a transition at an angle of approximately 50°, while no 
transition is predicted by Makkonen ’s results or in the experimental data. 
Similar transition locations and location of the maximum heat transfer co- 
efficient are predicted by both analytical methods for Reynold’s numbers 
of 2.8 x 10 5 and 8.8 x 10 s . However, the magnitude of the maximum heat 
transfer coefficient is over-predicted by the method used in LEWICE. 

Since the current method uses only a single roughness parameter, there 
is some question whether the maximum probable roughness or the equiva- 
lent sand grain roughness height should be used in comparisons with exper- 
imental data. A comparison of the analytical and experimental results with 
k a = 0.9 mm in LEWICE is shown in Figure B-5. These results compare 
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more favorably with the method of Makkonen and with experimental data 
than those in Figure B-4, but still over-predict the maximum value of the 
convective heat transfer coefficient. 

Measurements of convective heat transfer have also been made on simu- 
lated wooden ice accretion shapes roughened with grains of sand (Reference 
B-9). Figures B-6a, b, c, and d show the comparisons of the values cal- 
culated by LEWICE with experimental data for 2, 5, and 15 min glaze 
ice accretions. In these cases, the maximum probable roughness height of 
0.00033 m was used in the calculations. These results show that the pre- 
dicted values of h c generally compare well with experiment in the leading 
edge region where much of the accretion process will take place. 
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Figure B-l: Identification of the boundary layer terms used 
in the integral boundary layer method. 



Figure B-2: Identification of the computational surface being 
evaluated by the integral boundary layer method. 


181 



/KT NODES 

/N 


STAGNATION POINT 
(S - 0.0) 



Figure B-3: Analytical representation of a typical geometry 



10 20 30 40 50 60 

ANGLE. DEG 


Figure B-4: Comparison of caluclated and experimental corrective 
heat transfer characteristics, LEWICE k, = 0.00135m. 
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Figure B-5: Comparison of calculated and experimental corrective 
heat transfer characteristics, LEWICE k , = 0.0009m. 
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Figure B-6: Comparison of calculated and experimental heat 
transfer coefficients on simulated glaze and rime ice shapes. 
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b. Five min glaze ice shape. 
Figure B-6: continued 
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c. Fifteen min glaze ice shape. 
Figure B-6: continued 


186 






- 2.20 


1.76 -1.32 -0.BB -0. 14 0.00 0.4-4- 0.88 1.32 1.76 

s/c 


d. Fifteen min rime ice shape. 
Figure B-6: Concluded 
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APPENDIX C 

FLOW FIELD CORRECTIONS 

As discussed in the text, there are two types of flow field corrections 
that are applied in LEWICE to allow more accurate particle trajectory 
calculations. The first correction is made to avoid errors that occur in the 
flow field close to the surface because the geometry is represented by discrete 
points connected by line segments instead of a smooth curve. These are 
called discretization errors. The second correction is made to avoid errors 
caused by convoluted ice shapes. Both of these corrections are performed by 
producing an imaginary or pseudo-surface which is used in the calculations 
instead of the actual ice surface. The corrections are described in the 
following sections. 

Correction for Discretization Errors Near the Body 

The potential flow computer program used in this study has a relatively 
large discretization error very close to the body. As shown in Figure C-l, 
taken from Reference C-l, the longitudinal and vertical velocities around 
the leading edge of the Joukowski airfoil are very irregular close to the airfoil 
but become smoother as A£ is increased. (In Figure C-l, is equivalent to 
the DSHIFT parameter in LEWICE). The velocity is computed for different 
constant values of separation distance from the body, as illustrated in the 
insert. Note the large peaks or oscillations in the flow field velocity near 
the surface. 

Figure C-2 illustrates the trajectories of three different particles released 
at very small separation distances, y 0 , upstream of the body. Note that 
the upper particle turns near the nose and flows backward, reverses direc- 
tion again, and flows around the body. The initial intermediate positioned 
particle trajectory crosses the lower particle trajectory near the body. It 
approaches very close to the surface somewhat downstream of the nose, 
and then departs off into the free-stream without impinging on the body. 
Finally, the lower particle trajectory impinges on the body. These erratic 
trajectories are caused by the strong perturbations in the flow field near 
the body. 

To overcome the effect of the discretization error near the body, an 
artificial impingement surface is generated by the computer program. This 


188 


surface is achieved by displacing the surface of the body a small increment 
DSHIFT in the upstream x direction. To generate the pseudo-impingement 
surface, the x-coordinates are increased by the quantity DSHIFT times 
the cosine of the segment angle while the y-coordinates remain constant. 
The resulting pseudo-impingement surface is displaced outward into the 
freestream as shown in Figure C-3. The value of DSHIFT is input by the 
user in namelist TRAJl. Commonly used values between .2 to .6 percent 
of the chord length encompass the major area of uncertainty in the flow 
field. 

The pseudo-impingement surface is used only for particle impingement 
calculations and does not influence the potential flow calculations. The 
pseudo-impingement surface is discarded once the particle trajectory calcu- 
lations are completed and is not present when a new ice surface is formed. 

C.2 Correction for Multiple Calculated Stagnation Points 

Glaze ice accretions are often characterized by the formation of two 
horns, as shown in Figure C-4a. The calculation of the flow field around 
these typical glaze ice shapes produces undesirable results such as multiple 
calculated stagnation points (surface velocity = 0.0). An example of such a 
calculation is shown in Figure C-4b. When more than one stagnation point 
is calculated, errors occur in the particle trajectory and thermodynamic 
portions of the code, causing an abnormal termination. Since the Douglas 
potential flow solution is the most feasible method for use in this code, a 
method to obtain a sufficiently accurate flow field solution is required when 
this situation occurs. 

Before the development of computer codes capable of applying the 
Navier-Stokes equations to the calculation of flow characteristics in recir- 
culation zones, the potential flow codes available had a problem similar to 
that in this application. One solution was to replace the boundary of the 
separation/reattachment zone with a pseudo-surface, as shown in Figure 
C-5. A similar approach is used in this application to smooth the groove 
found at the stagnation point of many ice accretions. 

Determination of the Locations of Stagnation Points 

The location of the stagnation point(s) is determined in subroutine 
STAG by checking the non-dimensional edge velocities calculated by S24Y 
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(VT) for a change in sign. These velocities will be negative on the lower 
surface and positive on the upper surface. The velocity on each segment is 
checked, and the number of each segment at which the sign of VT changes 
is stored as a stagnation point. If more than one sign reversal is found, 
it is assumed that the flow field is inaccurate, and the user will be given 
the option of selecting one of the calculated stagnation points to use in 
the remaining calculations, or to form a pseudo-surface over the calculated 
stagnation points. 

The algorithm to form the pseudo-surface is located in subroutine 
PSURF. The computer prompts and responses used to form the surface 
can be found in Example 2 of Chapter 7. The purpose of the remainder of 
this appendix is to demonstrate the accuracy of this method by comparing 
experimental data to surface velocity profiles calculated using a pseudo- 
surface. The geometries used for these comparisons are simulated 2- and 
5-min glaze ice accretions and a 15-min rime ice accretion. These ice shapes 
are made of wood and simulate actual ice accretions formed on a cylinder. 
The surface velocity measurements were obtained from surface pressure 
data for each of the geometries. 

Recall that Figure C-4a shows a 5-min glaze ice accretion for which two 
stagnation points were calculated. The corresponding calculated surface 
velocity profile is shown in Figure C-4b. A pseudo-surface, shown in Fig- 
ure C-6a, was placed over the glaze accretion, and the resulting calculated 
surface velocity compared to experimental data. Note that the calculated 
velocity profile now has only one stagnation point, and that the calculated 
and experimental velocities compare very well near the stagnation point. 
Additional applications of the pseudo-surface technique are shown in Fig- 
ures C-6b and c. In all cases, the presence of the pseudo-surface produced 
smoother velocity profiles which more closely matched experimental results. 

This method, while capable of improving the accuracy of the flow field 
solution, is an approximation to the true solution and must be treated as 
such. For example, Figures C-7a and b show a 15-min glaze ice accretion 
and the corresponding calculated surface velocity. In this case, ten stag- 
nation points were calculated. A pseudo-surface, shown in Figure C-8a, 
was placed over the glaze ice accretion, and the calculations were repeated. 
As shown in Figure C-8b, the calculated velocity profile is much smoother 
than that for the actual ice shape, and the velocity in the groove is slightly 
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over-predicted. While this solution would allow the calculation of the ice 
accretion to continue, it illustrates that the pseudo-surface technique can 
produce a very rough approximation to the actual flow field. 
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C-l: Discretization error in the flow field near the nose of the Joukowski airfoil. 
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(3) FIVE NINUTE GLAZE ICE ACCRETION WITH GLAZE HORNS. 



(b) MULTIPLE CALCULATED STAGNATION POINTS FOR A 5 MINUTE GLAZE ICE SHAPE. 


C-4: Simulated five min glaze ice accretion on a cylinder and the corresponding 

calculated surface velocity profile. 
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a. Simulated five min glaze ice shape 

C-6: Comparison of experimental surface velocities on simulated ice shapes 
to those calculated for the ice shape with pseudo-surface. 
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b. Simulated two min glaze ice shape 
C-6: continued 
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c. Simulated fifteen min rime ice shape 
C-6: Concluded 
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b. Calculated surface velocity over the ice shape 
C-7: Concluded. 
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C-8a: Application of the pseudo-surface technique on a simulated 
fifteen min glaze ice accretion. 
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C-8b: Calculated surface velocity over ice shape with pseudo-surface. 
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APPENDIX D 

COMPLETE INPUT TO THE POTENTIAL FLOW CODE 
(S24Y) 

As noted in Section 5.1, the input to the potential flow code was simpli- 
fied for use in LEWICE. Many of the generalities in S24Y are not applicable 
to LEWICE, and, therefore, are set to default values in the program (SUB- 
ROUTINE SETUP). SUBROUTINE SETUP reads the simplified input 
from unit 35, and, after assigning the default values, writes to unit 45 in 
the proper S24Y input format. Since the original input to the code is still 
used, it is necessary to define all of variables used in this input file. The 
following is the description of the S24Y input file written to unit 45. 

D.l Potential Flow Code Input Cards 

Card 01 Program Control Card 

(3(11, 2X), IX, 7A4, 5X, 9(11, 2X), IX, II) 

Column Code Format Description 

01 ID Body Number 

04 ISV Flag to control the saving of the 

geometry for future use by the 2-D 
program 

= 0 Do not save data 
= 1 Save the input geometry data for 
future use 

07 ILIFT Lift Control Flag 

= 0 This is not a lifting body 
= 1 This is a lifting body 

11-38 TTITLE Body Description 
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IPARA 


IFIRST 


ISECND 


ITR 


INORM 


Element Geometry Flag 
= 0 Linear Elements 
= 1 Parabolic Elements 

First-order Terms Flag 
— 0 No first-order terms 
= 1 First derivative term 
= 2 Curvature term 
= 3 Both first-order terms 

Second-order terms flag 
= 0 No second-order terms 
= 1 Second derivative term 
= 2 Curvature squared term 
= 3 Both second-order terms 

Geometry transformation card 
= 0 Transformation card will not 
be input 

= 1 Geometry transformation card 
will be input 

= 2 Ellipse generation. Ellipse 
generation card will be input. 
Transformation card will not 
be input. 

= 3 Ellipse generation. Ellipse 
generation card will be input. 
Transformation card will be input. 

Geometry Normalization Card 
= 0 Geometry will not be normalized 
= 1 All of the geometry data 
(x and y) will be divided by 
the chord length before use by 
the potential flow program 
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IBOD 


Body Disposition Flag 
This flag together with the IDOLD 
parameter controls the sequence of 
the potential flow analysis part of 
the program. With the use of these 
two flags and the ISV parameter, it 
is possible to perform a variety of 
multi-element problems with a 
minimum of input data. For normal 
useage when all of the geometry data 
are input, only the IBOD = 1 and = 2 
are used. 

= 1 New geometry is being used. 

The storage of geometry data 
for the potential flow solution 
will start with this body. 

= 2 New geometry is being input, but 
this is not the first body. 

This body will be added to the 
sequence of body data already input. 
= 3 New geometry is being input, but 
it is to be added to an old 
sequence of data. 

= 4 All previously saved geometry will 
be used 

= 5 The geometry for this body will 
be selected from the previously 
saved data (body IDOLD will be 
selected). This selected body 
will be added to the current string 
of bodies. 

= 6 Previously saved geometry data 
will be used with the body 
number indicated by the IDOLD 
parameter removed from the solution 
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IDOLD 


65 


IPVOR 


Old Body ID Number 
This parameter is used in 
conjunction with the IBOD parameter 
in selecting which previously saved 
shape is to be retrieved as the 
present body 

Vorticity Distribution Flag 
= 0 Use constant vorticity between 
the body elements 

= 1 Use a variable vorticity distribution 
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LAST 


Last Body Flag 
= 0 This is not the last body. 

After this body is input, the 
program will return to read another 
Body Title and Control Card for the 
next body. 

= 1 This is the last body. 
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ITYPE 


Card Type Flag 
= 1 


X-Coordinate Cards (6F12.7, 2X, II, IX, II, IX, II) 


Column 

Code 

Format 

Description 

1-12 

X(l) 

6F12.7 

x-coordinates of the body geometry. 

13-24 

X(2) 


Up to six points may be on each card 

25-36 

X(3) 


depending upon how the INO flag is set. 

37-48 

X(4) 



49-60 

X(5) 



61-72 

X{6) 



75 

INO 

11 

Number of data points on the card 
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(Maximum of six) 

77 ISTAT II Last Card Flag 

= 0 This is not the last 

x-coordinate card. More cards will follow. 
= 1 This is the last x-coordinate card. 

79 ITYPE II Card Type Flag 

= 3 

Y-Coordinate Cards (6F12.7, 2X, II, IX, II, IX, II) 


Column 

Code 

Format 

Description 

1-12 

Y(I) 

6F12.7 

y-coordinates of the body geometry. 

13-24 

Y(2) 


Up to six points may be input on 

25-36 

Y(3) 


each card depending upon how the 

37-48 

Y(4) 


INO flag is set. 

49-60 

Y(5) 



61-72 

Y(6) 



75 

INO 

11 

Number of data points on the card 




(Maximum of six) 

77 

ISTAT 

11 

Last Card Flag 


= 0 This is not the last 

y-coordinate card. More cards 
will follow. 


=1 This is the last y-coordinate 
card 

79 ITYPE II Card Type Flag 

=4 
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Flow Title Card (15A4, 11X, II) 


Column 

Code 

Format 

Description 

1-60 

FTITLE 

15A4 

Title 

72 

ITYPE 

11 

Card Type Flag 
= 8 

Flow Control Card 

(11, 4X,F10.5,2X, 11, 2X,F10.5,5(4X.I1),9X, 11, 3X, 11,12,11) 

Column 

1 

Code 

INCLT 

Format 

11 

Description 
ci, a Flag 

= 0 Airfoil angle of attack, a, is 
input. 

6-15 

CLT 

F10.5 

= 1 Total lift coefficient, C t , is 

Value of the angle of attack or lift 
coefficient depending on how the 
INCLT flag was set 

18 

ICHORD 

11 

Reference Length Flag 
= 0 The reference length used to 
calculate C t is set = 1.0 
= 1 The reference length used to 

calculate C t will be input as CCL 

21-30 

CCL 

F10.5 

The input value for the reference 
length used in calculating C t 
(generally the airfoil chord) 

35 

IND 

11 

Individual Solution Flag 


S24Y is capable of calculating the 
potential flow about up to 6 bodies 
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and then superimpose the results of 
each. The possible values of IND are 
as follows: 

= 0 Edge velocities are not calculated 
for each body 

= 1 Edge velocities are calculated for 
each body 

In LEWICE, only one body is input and the 
edge velocities are always required. 
Therefore, IND = 1. 

Matrix Solution Method Control Flag 
= 0 Use routine SOLVIT for the 

matrix solution (used when a very 
large number of points have been input) 
= 1 Use routine QUASI for the matrix 
solution 

= 2 Use routine MISl for the matrix 
solution. The maximum number 
of geometry points = 101. If 
the number of points is greater 
than 101, the program will 
automatically shift to use SOLVIT. 

Off-body Calculation Flag 
= 0 Off-body points will not be 
calculated 

= 1 Off-body points will be calculated 

Non-uniform Flow Flag 
= 0 Non-uniform flow is not input 
= 0 Non-uniform flow will be input. 

The number of flows input = NONU 
(maxim um of six). When this option 
is used, the program sets ISOL = 1. 
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NBNU 

11 

The number of bodies for which 
the non-uniform flows are input 

65 

IPRINT 

11 

Print/punch Flag 
= 0 Normal output 
= 2 Print the individual matrices 
= 7 Punch the output on cards 

70 

MORE 

11 

Last Case Flag 


= 0 This is the last case 
= 1 This is not the last case 
Another set of Flow Title and 
Flow Control cards (and any 
non-uniform or off-body cards) 
will be expected after this case 
is completed 

70-71 IFILL 12 Parabolic Integration Flag 

S24Y calculates the forces and 
moments acting on the body using 
both a trapezoidal and parabolic 
integration of the calculated 
coefficient. The results of 
the trapezoidal calculations 
are always output. The value 
of IFILL determines whether the 
parabolic results are printed. 

72 ITYPE II Card Type Flag 

= 9 

Off-Body Title and Control Card 

(11, 9X, 7A4, 12X, 2(2X, II), 2(5X, II), 2X, 12) 
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Column Code 


Format 


Description 


1 ID II 


11-38 TITLE 7A4 

53 ITR II 


56 INORM II 


62 IDOLD II 


68 LAST II 


Identification number for this 
group of off-body points. Off-body 
points are read in groups of up to 
100. There is no limit on the number 
of groups. 

Title or description for this 
group of off-body points 

Coordinate Transformation Flag 

See the ITR parameter on the Body Title 

and Control Card 

Coordinate Normalization Flag 
= 0 Off-body coordinates will not 
be normalized 

= 1 Normalize the coordinates by 
the input chord or by the chord 
for the body with ID = IDOLD 

Body Selection Flag for Normalizing 
Off-body Points 

= 0 Use the input chord to normalize 
the off-body points 
= 0 Use the chord for body with 
ID = IDOLD to normalize the 
off-body points 

Off-body Group Termination Flag 
= 0 Another group of off-body 
points will follow this group 
= 1 This is the last group of off- 
body points 
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71-72 


ITYPE 12 


Card Type Flag 
= 21 


Off-Body X-Coordinate Cards (6F12.7, 2X, II, IX, II, IX, II) 


Column 

Code 

Format 

Description 

1-12 

13-24 

25-36 

37-48 

61-72 

X(l) 

X(2) 

X(3) 

X(4) 

X(6) 

6F12.7 

x-coordinates of the off-body 
points. Up to six points may be 
input on each card depending upon 
how the INO flag is set. 

75 

INO 

11 

Number of data points on the card 

77 

ISTAT 

11 

Last Card Flag 
= 0 This is not the last 


x-coordinate card. More cards 
will follow. 

= 1 This is the last x-coordinate 
card 

79 ITYPE II Card Type Flag 

= 3 

Off Body Y-Coordinate Cards (6F12.7, 2X, II, IX, II, IX, II) 


Column 

Code 

Format 

Description 

1-12 

Y(l) 

6F12.7 

y-coordinates of the off-body 

13-24 

Y(2) 


points. Up to six points may be 

25-36 

Y(3) 


input on each card depending upon 

37-48 

Y(4) 


how the INO flag is set. 

49-60 

Y(5) 



61-72 

Y(6) 


• 

75 

INO 

11 

Number of data points on the card 
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(Maximum of six) 


77 ISTAT II Last Card Flag 

= 0 This is not the last y-coordinate 
card. More cards will follow. 

= 1 This is the last y-coordinate 
card 

79 ITYPE II Card Type Flag 

= 4 


D.2 Sample Input 

The input to the S24Y code using the formats described in the previous 
section is written to unit 45 in subroutine SETUP. An example of the 
S24Y input file on unit 45 is shown in Figure D-l. For normal operation of 
LEWICE, this file should not have to be examined or modified. 
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I 0 1 NACA 0012 


1 3 

3 0 0 1 

oil 1 



1.0000000 

0.9899999 

0.9800000 

0.9700000 

0.9500000 

0.9250000 

6 

0 3 

0.9000000 

0.8750000 

0.8500000 

0.8250000 

0.8000000 

0.7750000 

6 

0 3 

0.7500000 

0.7250000 

0.7000000 

0.6750000 

0.6500000 

0.6250000 

6 

0 3 

0 .6000000 

0.5750000 

0.5500000 

0.5250000 

0.5000000 

0.4750000 

6 

0 3 

0 .4500000 

0.4250000 

0.4000000 

0.3750000 

0.3500000 

0.3250000 

6 

0 3 

0.3000000 

0.2750000 

0.2500000 

0.2250000 

0.2000000 

0.1750000 

6 

0 3 

0.1500000 

0.1250000 

0.1000000 

0.0900000 

0.0800000 

0.0700000 

6 

0 3 

0.0600000 

0.0500000 

0.0450000 

0.0400000 

0.0350000 

0.0300000 

6 

0 3 

0.0250000 

0.0200000 

0.0150000 

0.0100000 

0.0075000 

0.0050000 

6 

0 3 

0.0037500 

0.0025000 

0.0022500 

0.0020000 

0.0017500 

0.0015000 

6 

0 5 

0.0012500 

0.0010000 

0.0008750 

0.0007500 

0.0006250 

0.0005000 

6 

0 3 

0.0003750 

0.0002500 

0.0001250 

0.0000000 

0.0001250 

0.0002500 

6 

0 3 

0.0003750 

0.0005000 

0.0006250 

0.0007500 

0.0008750 

0.0010000 

6 

0 3 

0.0012500 

0.0015000 

0.0017500 

0.0020000 

0.0022500 

0.0025000 

6 

0 3 

0.0037500 

0.0050000 

0.0075000 

0.0100000 

0.0150000 

0.0200000 

6 

0 3 

0.0250000 

0.0300000 

0.0350000 

0.0400000 

0.0450000 

0.0500000 

6 

0 3 

0.0600000 

0.0700000 

0.0800000 

0.0900000 

0.1000000 

0.1250000 

6 

0 3 

0.1500000 

0.1750000 

0.2000000 

0.2250000 

0.2500000 

0.2750000 

6 

0 3 

0.3000000 

0.3250000 

0.3500000 

0.3750000 

0.4000000 

0.4250000 

6 

0 3 

0 . 4500000 

0.4750000 

0.5000000 

0.5250000 

0.5500000 

0.5750000 

6 

0 3 

0.6000000 

0.6250000 

0.6500000 

0.6750000 

0.7000000 

0.7250000 

6 

0 3 

0.7500000 

0.7750000 

0.8000000 

0.8250000 

0.8500000 

0.8750000 

6 

0 3 

0.9000000 

0.9250000 

0.9500000 

0.9700000 

0.9800000 

0.9899999 

6 

0 3 

1.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

1 

1 3 

0.0000000 

- 0.0018740 

- 0.0036110 

- 0.0052390 

- 0.0082510 

- 0.0117000 

6 

0 4 

- 0.0149080 

- 0.0179530 

- 0.0208880 

- 0.0237390 

- 0.0265150 

- 0.0292210 

6 

0 4 

- 0.0318550 

- 0.0344110 

- 0.0368870 

- 0.0392810 

- 0.0415850 

- 0 . 0437950 

6 

0 4 

- 0 . 0459040 

- 0.0479060 

- 0.0497930 

- 0.0515600 

- 0.0531980 

- 0.0546980 

6 

0 4 

- 0.0560510 

- 0.0572430 

- 0.0582620 

- 0.0590900 

- 0.0597070 

- 0.0600940 

6 

0 4 

- 0.0602260 

- 0.0600760 

- 0.0596120 

- 0.0587940 

- 0.0575700 

- 0.0558790 

6 

0 4 

“ 0.0536360 

- 0.0507320 

- 0.0470040 

- 0.0452280 

- 0.0432520 

- 0.0410380 

6 

0 4 

“ 0.0385350 

- 0.0356740 

- 0.0340820 

- 0.0323620 

- 0.0304960 

- 0.0284620 

6 

0 4 

- 0.0262300 

- 0.0237260 

- 0.0207970 

- 0.0172480 

- 0.0150780 

- 0 . 0123860 

6 

0 4 

- 0 . 0106990 

- 0.0086200 

- 0.0081360 

- 0.0076230 

- 0.0070750 

- 0.0064850 

6 

0 4 

- 0.0058470 

- 0.0051460 

- 0.0047660 

- 0.0043630 

- 0.0039310 

- 0.0034620 

6 

0 4 

- 0.0029450 

- 0.0023560 

- 0.0016320 

0.0000000 

0.0016320 

0 . 0023560 

6 

0 4 

0.0029450 

0.0034620 

0.0039310 

0.0043630 

0.0047660 

0.0051460 

6 

0 4 

0.0058470 

0.0064850 

0.0070750 

0.0076230 

0.0081360 

0.0086200 

6 

0 4 

0.0106990 

0.0123860 

0.0150780 

0.0172480 

0.0207970 

0.0237260 

6 

0 4 

0.0262300 

0.0284620 

0.0304960 

0.0323620 

0.0340820 

0.0356740 

6 

0 4 

0.0385350 

0.0410380 

0.0432520 

0.0452280 

0.0470040 

0.0507320 

6 

0 4 

0.0536360 

0.0558790 

0.0575700 

0.0587940 

0.0596120 

0.0600760 

6 

0 4 

0.0602260 

0.0600940 

0.0597070 

0.0590900 

0.0582620 

0.0572430 

6 

0 4 

0.0560510 

0.0546980 

0.0531980 

0.0515600 

0.0497930 

0.0479060 

6 

0 4 

0.0459040 

0.0437950 

0.0415850 

0.0392810 

0.0368870 

0.0344110 

6 

0 4 

0.0318550 

0.0292210 

0.0265150 

0.0237390 

0.0208880 

0.0179530 

6 

0 4 

0.0149080 

0.0117000 

0.0082510 

0.0052390 

0.0036110 

0.0018740 

6 

0 4 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

1 

1 4 

iCA 0012 





8 



4.00000 0 0.00000 1 

0 10 

0 

0 0 19 



NACA 

0012 



0 0 0 

1 1 



- 0.0200000 

- 0.0200000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

2 

1 3 

1 . 0000000 

0.5000000 

0.0000000 

0.0000000 

0.0000000 

0.0000000 

2 

1 4 


D-l : Input to S24y written to unit 45. 
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APPENDIX E 

CALCULATION OF THE LOCAL EDGE VELOCITIES 

Another pseudo-surface, related to the one created for particle impinge- 
ment calculations, is created to determine non-dimensional edge velocities. 
This pseudo-surface is formed by moving each point on the body a dis- 
tance DSHIFT along the outward normal to the surface. The resulting 
pseudo-surface is shown in Figure E-l. The non-dimensional edge velocities 
that were originally computed in subroutine FLOWS on the “true” surface 
are recalculated in subroutine READIN on the pseudo-surface. (This re- 
calculation step can be omitted with only a small change to subroutine 
READIN.) The purpose of creating a pseudo-surface is to shift the surface 
outward into a region where the velocity profile is smoother, thus avoiding 
numerical problems that may occur when velocities must be computed on 
an irregular ice surface. This method will produce higher stagnation ve- 
locities but will affect the ice shape only in a relatively small region of the 
airfoil unless a value of DSHIFT larger than recommended in this manual 
is used. 

The non-dimensional edge velocities that are recalculated on the pseudo- 
surface are the ones used in the evaluation of the thermodynamic charac- 
teristics of the ice surface. These non-dimensional velocities, V e ', are given 
by the equation 


V' e = £- {E-l) 

The values of V' are initially written to unit 30 in subroutine READIN 
and then read in subroutine VEDGE, where they are used to calculate the 
dimensional local edge velocities. The boundary layer edge velocities, V e , 
temperature, T e , and pressure, P e , are calculated from V' e using the isen- 
tropic equations with the Karman-Tsien compressibility correction. These 
parametes are assigned the following variable names in the computer code: 


Local static pressure, (P e ) PE(I) 

Local static temperature, ( T t ) TE(I) 
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Local dimensional velocity, {V t ) VE(I) (compressible) 

Local non-dimensional velocity, (V,') VT(I) (incompressible) 


Each of the parameters are arrays with (I) denoting the segment number. 
This appendix describes the equations used to calculate these flow proper- 
ties. 

The free stream mach number, total temperature, and total pressure 
are first calculated using the following isentropic equations: 


M - V °° 

°° 20.05v^7 

(■ E-i ) 

M 2 

T t = T.( 1 + £*) 

(E — 3) 

1/2 

P r = P.(l + ^) 3 - 5 

(E- 4) 

An incompressible pressure coefficient is calculated 
dimensional edge velocities using the following equation: 

from the non- 

= 1 - W 

(E- 5) 


This incompressible pressure distribution around the body is then corrected 
for compressibility using the Karman-Tsien compressibility correction given 
by the following equation: 


c 


Pc 


(1 - MlY” + 


C P.n« 

2 


1 . + (1 - 


(E-6) 


The local pressure is calculated using this corrected pressure coefficient and 
the following equation: 
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P, = F.(l + |c,.Mi) 


(£-7) 


where P t is the free stream, static pressure. If a local static pressure is 
calculated to be greater than the total pressure, the local static pressure is 
set equal to the total pressure. 

The local mach number is calculated using the following isentropic re- 
lationship between the local static and total pressures: 


M e = 



[E — 8) 


The local static temperature is determined using the following isentropic 
relation: 


r. = r r ^i + ^ (£-9) 

The dimensional local edge velocity is then calculated using the isentropic 
equation for the speed of sound in air and the local Mach number, e.g., 

V e = M e (20.05y/T t ) ( E - 10) 


The values of VE, TE, and PE for each body segment are then written to 
unit 30. 


% 
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E-l : Illustration of the pseudo-surface 


APPENDIX F 

EMPIRICAL RELATIONSHIP FOR THE EQUIVALENT 
SAND GRAIN ROUGHNESS HEIGHT 

As discussed in Section 5.2, the size, shape, and type of ice accretion 
that is formed is dependent upon the convective heat transfer rate from the 
ice surface. The integral boundary layer method, described in Appendix B, 
requires that a surface roughness height be specified to identify transition 
to turbulent flow and evaluate the convective heat transfer characteristics 
of the rough surface. 

Increasing the size of the roughness elements will increase the calculated 
convective heat transfer coefficient. Therefore, a calculated ice shape can 
change quite drastically when various values of roughness element height are 
specified. Figure F-l shows the effect of varying the input value of the sand 
grain roughness height, k„ on the calculated ice accretion shape. When 
the value of k, is smaller, the amount of heat removed from the surface by 
convection is reduced. Therefore, the surface temperature is not lowered 
sufficiently to allow all of the ice to freeze on impact and a glaze accretion is 
formed. Incrementally increasing k, increases the convective heat transfer, 
thereby increasing the fraction of impinging water that freezes on contact. 
As a result, the accretions begin to take on more characteristics of rime 
ice. At some value of k„ sufficient heat will be removed to freeze all of the 
water on impact, and a complete rime ice accretion is formed. 

Empirical correlations that can be used to evaluate the effect of rough- 
ness exist in the literature (Reference F-l and F-2). Unfortunately, these 
correlations are usually applicable only for a specific, well-defined type of 
roughness element. None of these correlations are directly applicable to the 
irregular surface roughness elements found on typical ice accretions. Fur- 
thermore, the size and shape of the roughness elements on an ice surface are 
dependent on the conditions at which the ice was formed. They can vary 
with the icing condition, surface location, and, since the ice shape changes 
with time, the icing time. While the effect of each of these parameters on 
the formation of the surface roughness elements is very complex, all analyt- 
ical ice accretion prediction methods must address the problem because of 
the strong influence on the calculated ice accretion shape shown in Figure 
F-l. 
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It has become standard practice in current ice accretion prediction 
methods to develop an empirical correlation for either the surface element 
roughness height (Reference F-3) or the convective heat transfer coefficient 
(Reference F-4). These correlations are developed by first predicting the ice 
shapes for a set of experimental ice shapes by changing the convective heat 
transfer coefficient (or roughness element height) to determine the value 
that yields the best agreement with experiment. Unfortunately, the corre- 
lation will depend on the computational algorithm and may not be directly 
applicable to any other analytical ice accretion prediction method. Since 
the timestepping procedure applied in LEWICE is unique among analyti- 
cal ice accretion prediction methods, an empirical correlation relating the 
surface roughness height to the icing condition had to be developed. 

The experimental data used to develop this correlation was obtained 
by Gent (Reference F-5). This data set, shown in Figures F-2, -3, and -4, 
show the effect of velocity, LWC, and static temperature on the shape of 
the ice accretion formed. Also shown in the figures are the ice accretion 
shapes, calculated by LEWICE, that best compare with the experimental 
shape and the corresponding value of k t . The accretion shown in Figure F- 
3a for LWC =0.5 g/m a was used as a baseline condition for the correlation. 
Therefore, all values of k t were divided by the value of k , for this case and 
plotted as a function of either velocity, LWC, or static temperature. The 
resulting data points, normalized by the airfoil chord of 0.3m are shown in 
Figures F-5, -6, and -7, respectively. 

A correlation relating the sand grain roughness height to the icing pa- 
rameter was formed by fitting the data on each of these plots with a least 
squares linear or quadratic curve fit. The curve calculated from each cor- 
relation is also shown in Figures F-5, -6, and -7. The equations are as 
follows: 

Velocity 

, k *( C -- = 0.4286 + 0.0044139(V oo ) (F - 1) 

k$j CJbate 


Liquid Water Content 

= 0.5714 + 0.2457(LWC) + 1.2571 (LWC) 3 (F - 2) 

K./C)ba»e 
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Static Temperature 


ki/c)base 


46.8384 



11.2037 


(F — 3) 


In all of these equations, the baseline value, is 0.00117. The ve- 

locity, LWC, and static temperature are input into the equations in m/sec, 
g/m* , and K, respectively. By relating them to a baseline value, the corre- 
lations calculate a multiplying factor accounting for the effect of velocity, 
LWC, and static temperature. The value of sand-grain roughness height to 
be input is calculated using the following equation: 


k. 


k 8 j c 


kjc 


kjc 

k 8 j c)b a$e 

Voo 

k& l Cubase m 

LWC 

k 9 / c)t 0#e 


k t j c)bate ® 


(F — 4) 


The use of these equations is best illustrated by a numerical example. 
Suppose that the ice shape for the following icing condition is to be deter- 
mined. 


Airfoil Chord 
Velocity 

Static Temperature 
LWC 


= 0.3m 

= 75.0 m/s 

= 255.0 K 9 

= 0.5 g/m 


The multiplying factors calculated from Equations (F-l), (F-2), and (F- 
3) are 0.7596, 1.0085, and 0.7401, respectively. Substituting these values 
into Equation (F-4), the final factor to be multiplied by the baseline value 
and the airfoil chord is 0.5670. Therefore, for a baseline value of 0.00117 
and chord of 0.3 m, the sand-grain roughness to be input into the code is 
0.000199 m. 
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As previously noted, these three correlations were developed to account 
for the effect of velocity, LWC, and static temperature on the surface rough- 
ness elements formed on an ice accretion. These three icing parameters were 
selected because of their obvious influence on the type of ice formed, and 
because a complete and consistent set of experimental data existed from 
which to form a correlation. The effects of droplet diameter, body geometry, 
static pressure, etc. have not been included. Also inherent in the baseline 
value of sand-grain roughness height are any characteristics unique to the 
facility in which the experimental ice accretion shapes were formed. These 
can include levels of free-stream turbulence, LWC and droplet diameter cal- 
ibrations, and possibly even the droplet size distribution produced by the 
spray nozzle. Comparisons with experimental ice accretion shapes have in- 
dicated that the baseline value of fc,/c)t a#e = 0.00117 yields good results for 
the facility of Reference F-l and the NASA Lewis Icing Research Tunnel 
(IRT) but would not be expected to be appropriate for all icing facilities 
and applications. A value for flight test data has not been determined but 
is expected to be less than 0.00117, primarily because of the lower level of 
free-stream turbulence encountered in flight. 

The surface roughness, while known to be a function of icing time and 
surface location, is currently specified to be constant for the entire icing 
time. Sufficient experimental data does not exist to produce meaningful 
correlations relating these parameters to the roughness element height. 
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RIME FEATHERS 






VELOCITY (Hz'S) 93.08 
TEMPERATURE <C> -19.37 
PRESSURE (KPA) 93.61 
HUMIDITY < X ) 100.00 
LHC (G/M--3) 0.96 
DROP DIAM (MICRONS) 37.50 
TIME (SECT 225.00 


F-l: Effect of varying the equivalent sand-grain roughness 
on the ice accretion calculated for a mixed icing condition 
on a NACA 0012 airfoil section at 0.0 deg angle of attack. 
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EXPERIMENTAL 


THEORETICAL 


A 



LHC - 0.25 9 /H 3 




LHC - 0.75 g/n 3 



LHC * 1.00 g/M 5 


VELOCITY < M/S ) 

1 27.00 

TEMPERATURE <C> 

- 12.60 

PRESSURE (KPA) 

90. ?6 

HUMIDITY <X> 

100.00 

DROP 0 1 AH (HICRONS) 

20.00 

TIME < SEC) 

120.00 


a. Angle of attack = 0.0° 

F-3: Effect of LWC on ice accretion shape. 
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b. Angle of attack = 4.0° 
F-3: Concluded. 
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F-4: Effect of static temperature on ice accretion shape 


229 




i.5 r 



0 50 100 150 200 

VELOCITY, M/S 

F-5: Emperical relationship for equivalent sand-grain roughness 
as a function of velocity. 


i 
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F-6: Emperical relationship for equivalent sand-grain roughness 

as a function of LWC. 
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F-7: Emperical relationship for equivalent sand-grain roughness 
as a function of static temperature. 
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